Ecology

Field experiment
Estimation of turtle catch per unit effort
We conducted our field experiment in February–April 2018 at two wetland complexes near Murray Bridge, South Australia, selecting two study sites at each complex (Supplementary Fig. S3). At each site, we estimated turtle population density using catch-per-unit-effort (CPUE; Supplementary Table S1). We conducted three 3-day rounds of turtle trapping using a combination of fyke and cathedral traps, baited with offal. Up to eight traps were deployed at a time. We calculated turtle CPUE by dividing the total number of turtles caught (regardless of species) by the total trap-hours. The number of trap-hours was similar across all four sites (average 1685 ± 7.6 SE total trap-hours).
Carp carcass decomposition
After the first and the second trapping rounds, we deployed whole carp carcasses at each site to measure carp decomposition rates depending on turtle accessibility. We placed each carp in a pre-weighed plastic box (340 × 230 × 120 mm), securing it with cable ties. Carp were made non-accessible to turtles in half of the deployments by covering the plastic boxes with 25 × 25 mm mesh (Supplementary Fig. S4). The mesh prevented turtle access to the carp, but was large enough to allow scavenging by crayfish (Cherax destructor) and other freshwater invertebrates. We tied each box to a brick and submerged the boxes around the four study sites ≥ 30 m away from each other, sunk at an average depth of 436 mm (± 13 SE). We used a total of 38 accessible and 40 non-accessible carp, split between our four study sites over two rounds (Supplementary Table S9). Every day, starting from day 2, the box and carp were weighed together with a digital scale. In all measurements, we calculated the wet mass of the carp by subtracting the box weight from the total weight. Carp carcasses were left in the wetlands for up to 10 days, or until they were fully consumed. All work was performed in accordance with DEWNR Permit M26663-1, PIRSA permits MP0085 and ME9902980, and The University of Sydney Animal Ethics Committee approval (project number 2017/1208), observing all relevant guidelines and regulations.
Statistical analysis
We analysed our data using RStudio 1.1.45633 (packages: “lme4” 1.1-2134, “MuMIn” 1.42.135). To assess whether turtles were important scavengers of our carcasses, we computed a linear mixed model testing whether turtle CPUE and carp access (yes/no) affected the rate of mass loss of the carp carcasses. The turtle CPUE values used were the average CPUE in the trapping round before and after each carp was deployed. We used the rate of mass loss per day as a dependent variable. We included the carp mass before deployment as an independent variable to account for initial mass variation, and we included study site as a random variable. We log-transformed all data before analysis. We assessed model fit by examining predicted versus residual and Q–Q plots, and testing the normality of residuals.
Mesocosm experiment
Turtle trapping and experiment procedure
We caught 20 adult male E. macquarii with fyke nets baited with offal at Hawkesview Lagoon, Albury, NSW, in November 2018. The E. macquarii captured at this site belong to the same genetic population as the E. macquarii trapped in South Australia36, therefore we expect behaviours to be similar between the two populations. We focussed on E. macquarii as this is the most common species in the Murray–Darling Basin, and fish carrion is an important part of its diet19,37. The turtles were transported by car to the Experimental Wetlands facility at Western Sydney University, in Richmond, NSW (Supplementary Fig. S5). This facility is comprised of 10 circular mesocosms (0.42 m depth × 2.1 m diameter) filled with 1,450 L of tap water. Each is an independent flow-through system where the water flow is regulated, and was maintained at 1998.6 ml/min (± 149.5 SE) throughout the experiment. Each mesocosm had two cement blocks for the turtles to bask on, and two plastic tunnels for shelter. The experiment was conducted for 40 days, therefore it is a short-term study (Supplementary Fig. S6). Upon arrival at the facility, we placed four adult male E. macquarii turtles in each of five random mesocosms, which means the experimental replication was 5. The remaining five mesocosms were controls and had no turtles. The four turtles comprised an average 5,376.6 g total biomass per pond, each being 3.46 m2. This would result in a biomass of 11,560 individuals/ha or 15,537 kg/ha on average. Kinosternon integrum has been estimated reaching densities of 20,000 individuals/ha in Sonora, Mexico, while Podocnemis vogli may reach 10,300 individuals/ha or 15,450 kg/ha in Venezuela, likely in temporary aggregations38. Emydura macquarii tend to congregate around food sources, therefore we considered four turtles per carp carrion as a realistic density. After 7 days of acclimation, we introduced one carp carcass to all mesocosms, and a second 6 days later. We used one ~ 1 kg carp at a time to simulate a density close to 3,144 kg carp/ha31. The turtles had continual access to the carp, which was their main food source throughout the experiment. The day all carp carcasses were fully eaten in all turtle mesocosms, we removed turtles from their mesocosms and released them at the point of capture. On the same day (day 10), we ended the data collection in their mesocosms, because any further change in water quality here would not have been related to carp decomposition. We continued the daily water quality measurements in the five control mesocosms until all carp were fully decomposed (day 32). This experimental design allowed us to collect water quality data without the need to add turtle food to the mesocosms, which would have biased our measurements once carp were removed from the turtle mesocosms.
We measured water temperature, dissolved oxygen, conductivity, turbidity, phosphate, and ammonia concentration in all mesocosms every morning from the day before the first carp introduction (see Supplementary Materials for equipment used). We also photographed the carcasses daily to estimate their decomposition rate based on a scale (Supplementary Table S10) designed after the decomposition stages described by Benninger et al.39 Due to the short transit time of fish matter in E. macquarii’s gut37, the effects of the turtles’ metabolic wastes on water quality are included in our experiment for carp 1. All work was performed in accordance with OEH Permit SL100401, DPI permit P09/0070-3.0, and Western Sydney University Animal Ethics Committee Animal Research Authority approval A12390, observing all relevant guidelines and regulations.
Statistical analysis
To assess whether the presence of turtles affected the decomposition of carp we computed a mixed linear model using the repeated measures PROC MIXED procedure using SAS (3.8 University Edition, SAS Institute Inc., Cary, NC, USA). For this model, days to total decomposition/removal was a dependent variable, turtle presence/absence was a fixed effect, and carp number (first or second) was a repeated fixed effect.
We used DO, conductivity, turbidity, phosphate, and ammonia to carry out a principal component analysis (PCA) using PROC PRINCOMP. We conducted a PCA because the parameters are a multivariate response and have potential to covary with each other, which would not be detected in univariate analyses. We considered a parameter loaded onto a PC when the absolute value of its eigenvector was > 0.300. If the same parameter loaded onto more than one PC, we considered it only on the PC where its eigenvector had a higher absolute value.
To test the effect of turtles scavenging on water quality, we computed general linear mixed models (GLMMs), using PCs as response variables, in PROC MIXED. We used a PC as response variable if its eigenvalue was greater than one (Kaiser criterion40). For each of these GLMMs, we included turtle presence (yes/no) and day number after the first carp introduction as fixed effects, water temperature and flow as covariates, and mesocosm ID as a random effect. We computed a model with full interactions first, and then, in absence of four- or three-way interactions, simplified the model to focus on main effects.
Finally, to test the effect of turtle scavenging on each water quality parameter, a GLMM was computed for each (logged) parameter that loaded onto a PC with eigenvalue > 1, i.e. dissolved oxygen, ammonia, turbidity, conductivity, phosphate. For these GLMMs, turtle presence and day number were fixed effects, water temperature and flow were covariates, and mesocosm ID was a random effect. More

Animals and housing conditions
The present study was conducted with 24 male C57BL/6J mice, purchased from a professional breeder (Charles River Laboratories, Research Models and Services, Germany GmbH, Sulzfeld, Germany) at the age of five weeks. Upon arrival, mice were housed in same-sex groups of 3 individuals per cage (Makrolon cages type III, 38 × 23 × 15 cm3), since in sub-adult male mice, the occurrence of escalated aggression is very unlikely. However, with the males becoming adult, the probability of escalated agonistic encounters increases. Therefore, at the age of nine weeks, mice were transferred to single housing conditions to avoid any escalated aggressive interactions. Please note that the question whether to house male laboratory mice singly or in groups is under ongoing discussion and there is still no “gold standard” regarding its solution. For current discussions about recommendations for male mouse housing see37,38. Cages were equipped with wood chips as bedding material (TierWohl Super, J. Rettenmaier & Söhne GmbH + Co.KG, Rosenberg, Germany), a wooden stick, a semi-transparent red plastic house (11.1 × 11.1 × 5.5 cm3, Tecniplast Deutschland GmbH, Hohenpeißenberg, Germany), and a paper tissue. Housing rooms were maintained at a reversed 12 h dark/light cycle with lights off at 8 a.m., a temperature of approximately 23 °C, and a relative humidity of about 50%. The animals had ad libitum access to water and food (Altromin 1324, Altromin Spezialfutter GmbH & Co. KG, Lage, Germany) until the beginning of the touchscreen training phase. From then on they were mildly food restricted to 90–95% of their ad libitum feeding weights in order to enhance their motivation to work for food rewards. As neither distinct negative effects of such a restricted feeding protocol39, nor an interference with judgement bias assessment17,18 could be detected in previous studies, we considered this method to not affect the emotional state of the mice itself. Weights were monitored on a daily basis using a digital scale (weighing capacity: 150 g, resolution: 0.1 g; CM 150-1 N, Kern, Ballingen, Germany).
In addition to the experimental animals, 16 group-housed adult female C57BL/6J mice and 5 single-housed adult male NMRI mice, purchased from Charles River Laboratories, were used to provide the test animals with social experiences.
Ethics statement
All procedures complied with the regulations covering animal experimentation within Germany (Animal Welfare Act), the EU (European Communities Council DIRECTIVE 2010/63/EU), and the fundamental principles of the Basel Declaration, and were approved by the local (Gesundheits- und Veterinäramt Münster, Nordrhein-Westfalen) and federal authorities (Landesamt für Natur, Umwelt und Verbraucherschutz Nordrhein-Westfalen “LANUV NRW”, reference number 84-02.04.2015.A441).
Experimental design
In this study, the effects of different social experiences on important correlates of animal emotions, comprising cognitive (judgement bias), behavioural (anxiety-like and exploratory behaviour) as well as physiological (stress hormone levels) measures, were investigated. The experiment comprised six phases: a handling phase, a training phase, a first cognitive judgement bias (CJB) test phase, an experience phase, a second CJB test phase, and a behavioural test phase (Fig. 1).
Figure 1

Experimental design. Mice were habituated to cup handling before they underwent daily training sessions until successful acquisition of the discrimination task. Afterwards, they were tested in the cognitive judgement bias (CJB) test. During the following phase, mice repeatedly made one out of three different experiences: mildly “adverse”, “beneficial”, or “neutral”. They were then tested for their CJB again. During this second test phase, a so-called reminder was presented before each test session with the aim to re-evoke the affective state the animals experienced during the treatment phase. On the last day of each CJB test phase, faecal corticosterone metabolite concentrations (FCMs) were assessed. Subsequently, animals were tested for anxiety-like behaviour. Again, they were presented with reminders before each behavioural test.

Full size image

During the handling phase starting at PND 69, mice were first habituated to cup handling for 5 days and thereafter underwent daily training sessions to learn the discrimination task required for CJB testing, starting at PND 76. Afterwards, the animals’ initial CJB was assessed (start test phase 1: PND 223 ± 77; for details on CJB training and testing see following section).
During a subsequent experience phase starting at PND 230 ± 77, mice were exposed to one of three different experiences, each comprising three group-specific encounters, classified as either mildly “adverse”, “beneficial”, or “neutral”. Encounters took place under red light between 2:45 p.m. and 4:35 p.m. on 3 different days, always separated by a gap day. The mildly “adverse experience” group (AE group, n = 8) repeatedly encountered a dominant opponent of the aggressive NMRI strain40, with each confrontation lasting maximally 10 minutes30,41. Confrontations were terminated in cases of high aggression. The “beneficial experience” group (BE group, n = 8) was repeatedly presented with freshly collected urine of an unfamiliar C57BL/6J female for 10 minutes31. To provide all subjects with comparable experiences, we controlled for the females’ oestrus state. Since the time of oestrus in mice is relatively short42, urine from non-oestrous females was used in order to keep the total number of involved females low. The “neutral experience” group (NE group, n = 8) served as a control group and was repeatedly placed into a novel cage containing clean bedding material for 10 min.
Following the experience phase, CJB was assessed again to investigate the influence of the respective experience on the animals’ judgement bias (start test phase 2: PND 237 ± 77). In this second test phase, a so-called reminder was presented immediately before each test session. These reminders were introduced to acutely re-evoke the affective state the mice experienced during the encounters of the treatment phase. Reminders took place immediately before each test session of the second CJB test phase. For this purpose, mice were placed into a cage (Makrolon type II cage; 22 × 16 × 14 cm3) filled with bedding for 3 min. For AE mice, another 25 ml of soiled bedding from the home cage of the last NMRI male encountered were added. For BE mice, the same was done with soiled bedding from the home cage of the last female of which urine had been presented.
On the last day of each CJB test phase, faeces samples were obtained to assess corticosterone metabolite (FCM) concentrations. Finally, animals underwent a battery of standard behavioural tests for anxiety-like behaviour and exploratory locomotion (elevated plus maze test (EPM), dark-light test (DL), and open field test (OF); start: PND 245 ± 77). Before each test session, a reminder was presented again.
The allocation of mice to the treatment groups was pseudo-randomised, so that balanced numbers of mice with different learning speeds were present in each group. The testing order of mice was randomised once before the first CJB test and subsequently maintained for the following CJB and behavioural test sessions. As reminders were provided immediately before CJB testing as well as before the subsequent behavioural tests, blinding of the experimenter was not possible.
The touchscreen-based cognitive judgement bias test
Procedure
The same apparatus as described previously was used28,36 (Bussey-Saksida Mouse Touch Screen Chambers, Model 80614, Campden Instruments Ltd., Loughborough, Leics., UK). Mice underwent daily touchscreen sessions at intervals of approximately 24 h on maximally 6 consecutive days. Before each session, each mouse was taken out of its home cage and weighed. In a red semi-transparent box (21 × 21 × 15 cm3) the animal was then transported to a separate room where it was placed into the touchscreen chamber. The session was started and ended after a maximum number of trials had been performed or after a training step-specific duration. All touchscreen sessions were conducted during the dark phase between 8.15 a.m. and 1 p.m.
Paradigm
The paradigm applied here was the same as described previously with minor modifications36. Briefly, mice were trained to distinguish between a positive and a negative condition (Fig. 2). The positive condition was signalled by a bar at the bottom (5 cm below upper edge) of the cue-presentation field, the negative condition by a bar at the top (1 cm below upper edge). Mice had to touch either the left or right touch field in response to the cues. A correct touch in the positive condition led to the delivery of a large reward (12 μl of sweet condensed milk, diluted 1:4 in tap water, in the following “SCM”). An incorrect touch resulted in the delivery of a small reward (4 μl of SCM). In the negative condition, correct touches led to the delivery of a small reward (4 μl of SCM), while incorrect touches resulted in a mild “punishment” (5 s time out and houselight on). Mice had to learn to touch the high-rewarded side in the positive condition and the small-rewarded side in the negative condition. The small-rewarded touch field was the same in both conditions. The association between condition and correct touch side was the same for each individual but counterbalanced between mice. For a detailed description of the training procedure please see the supplementary material. After successful training, animals underwent CJB testing. The two cognitive bias test phases took place on five consecutive days each. During each CJB test session, three types of ambiguous cues, interspersed with the learned reference cues, were presented. These were bars at three intermediate positions: near positive (NP, 4 cm below upper edge), middle (M, 3 cm below upper edge) and near negative (NN, 2 cm below upper edge). Touches in response to these ambiguous cues resulted in a neutral outcome (neither a reward nor a “punishment”). The animals’ judgements made in response to these cues indicated whether they interpreted them according to the positive (“optimistic” response) or negative (“pessimistic” response) reference cue, serving as a measure of CJB.
Figure 2

Touchscreen-based cognitive bias paradigm. Mice were trained to distinguish between bars displayed at the top (negative condition) or bottom (positive condition) of a central field of a touchscreen. In this example, mice learned to touch right for a large reward during the positive condition and to touch left for a small reward during the negative condition (the association between positive/negative cue and the correct touch side was counterbalanced across mice). During the test, mice were presented with cues displayed at three intermediate positions: near positive, middle and near negative. The relative number of “optimistic”-like responses to these ambiguous conditions served as outcome measures of the animals’ cognitive judgement bias. Figure adopted from Krakenberg et al. (2019) with permission from Elsevier36.

Full size image

Each test session comprised 54 trials. Per session, each type of ambiguous cue was presented twice, interspersed with 48 training trials. Per test phase, each mouse was presented with each ambiguous cue ten times and each trained cue 120 times.
Behavioural measures
Responses to ambiguous cues served as a measure of the animals’ CJB. Touches according to the positive condition were defined as “optimistic” choices, touches according to the negative condition were defined as “pessimistic” choices. Out of all responses per condition, a “choice score” was calculated as previously28,36 according to the following formula:

$$ Choice,Score = frac{{N,choices ( {text{“}}optimistic{text{”}} ) – N,choices ({text{“}}pessimistic{text{”}})}}{ N,choices ({text{“}}optimistic{text{”}} + {text{“}}pessimistic{text{”}})} $$

The choice score could range between − 1 to + 1. Higher scores indicated a higher proportion of “optimistic” choices and consequently a relatively positive CJB compared to lower scores. Please note that choice scores are not an absolute, but a relative measure of CJB and that the term was chosen for the sake of intuitiveness.
Anxiety-like behaviour and exploratory locomotion
Mice were tested in three tests on anxiety-like behaviour and exploratory locomotion in the following order: the elevated plus-maze test (EPM), the dark-light test (DL) and the open field test (OF). The sequence of tests followed recommendations to schedule tests that are more sensitive to previous experience at the beginning of such a battery, and to conduct potentially more stressful tests towards the end43,44. Tests were carried out at intervals of at least 48 h and were performed in a room different from the housing room between 12:45 p.m. and 3:35 p.m. Test equipment was cleaned with 70% ethanol between subjects. Behaviour was recorded with a webcam (Logitech Webcam Pro 9000) and the animals’ movements during the EPM and OF were automatically analysed by the video tracking system ANY-maze (ANY-maze version 4.99, Stoelting Co., Wood Dale, IL, USA). Videos of the DL were analysed manually by an experienced observer (Sophie Siestrup). For apparatus descriptions and details about testing procedures see supplementary material.
Faecal corticosterone metabolites
The basal levels of adrenocortical activity of the subjects were monitored non-invasively by measuring faecal corticosterone metabolites45,46,47 (FCMs). Faeces samples of each individual were collected on the last day of the first CJB test week (= before the experience phase) and on the last day of the second CJB test week (= after the experience phase). During the dark phase, a peak of FCMs can be found in the faeces 4–6 h after the exposure to a stressor45. For this reason, faeces samples were collected 5.5–8.5 h after an individual finished CJB testing to ensure that faeces collection could be finished in the dark phase. For sample collection, mice were placed in Makrolon cages type III with a thin layer of bedding material and clean enrichment items as present in the home cage. Water was available ad libitum. After the sampling period of 3 h, mice were transferred to novel clean cages together with the enrichment items. All faeces produced during this time were collected and frozen at − 20 °C. Faecal samples were dried and homogenised and aliquots of 0.05 g were extracted with 1 ml of 80% methanol. Samples were then analysed using a 5α-pregnane-3β, 11β, 21-triol-20-one enzyme immunoassay (for details see45,46). Intra- and inter-assay coefficients of variation were More

1.
Graham, C. H., Moritz, C. & Williams, S. E. Habitat history improves prediction of biodiversity in rainforest fauna. Proc. Natl Acad. Sci. USA 103, 632–636 (2006).
ADS  CAS  PubMed  Google Scholar 
2.
Jansson, R. Global patterns in endemism explained by past climatic change. Proc. R. Soc. B 270, 583–590 (2003).
PubMed  Google Scholar 

3.
Rosenzweig, M. L. Species Diversity in Space and Time. (Cambridge University Press, Cambridge, 1995).

4.
Rosauer, D. F. & Jetz, W. Phylogenetic endemism in terrestrial mammals. Glob. Ecol. Biogeogr. 24, 168–179 (2015).
Google Scholar 

5.
Dynesius, M. & Jansson, R. Evolutionary consequences of changes in species’ geographical distributions driven by Milankovitch climate oscillations. Proc. Natl Acad. Sci. USA 97, 9115–9120 (2000).
ADS  CAS  PubMed  Google Scholar 

6.
Sandel, B. et al. The influence of late quaternary climate-change velocity on species endemism. Science 334, 660–664 (2011).
ADS  CAS  PubMed  Google Scholar 

7.
Nathan, R. et al. Spread of North American wind-dispersed trees in future environments. Ecol. Lett. 14, 211–219 (2011).
PubMed  Google Scholar 

8.
Malcolm, J. R., Markham, A., Neilson, R. P. & Garaci, M. Estimated migration rates under scenarios of global climate change. J. Biogeogr. 29, 835–849 (2002).
Google Scholar 

9.
Parmesan, C. & Yohe, G. A globally coherent fingerprint of climate change impacts across natural systems. Nature 421, 37–42 (2003).
ADS  CAS  PubMed  Google Scholar 

10.
Donoghue, M. J. A phylogenetic perspective on the distribution of plant diversity. Proc. Natl Acad. Sci. USA 105, 11549–11555 (2008).
ADS  CAS  PubMed  Google Scholar 

11.
McLachlan, J. S., Clark, J. S. & Manos, P. S. Molecular indicators of tree migration capacity under rapid climate change. Ecology 86, 2088–2098 (2005).
Google Scholar 

12.
Thuiller, W. et al. Consequences of climate change on the tree of life in Europe. Nature 470, 531–534 (2011).
ADS  CAS  PubMed  Google Scholar 

13.
Gonzalez-Orozco, C. E. et al. Phylogenetic approaches reveal biodiversity threats under climate change. Nat. Clim. Change 6, 1110–1114 (2016).
ADS  Google Scholar 

14.
Pollock, L. J., Thuiller, W. & Jetz, W. Large conservation gains possible for global biodiversity facets. Nature 546, 141–144 (2017).
ADS  CAS  PubMed  Google Scholar 

15.
Forest, F. et al. Preserving the evolutionary potential of floras in biodiversity hotspots. Nature 445, 757–760 (2007).
ADS  CAS  PubMed  Google Scholar 

16.
Weinstein, B. G. et al. Taxonomic, phylogenetic, and trait beta diversity in South American hummingbirds. Am. Nat. 184, 211–224 (2014).
PubMed  Google Scholar 

17.
IPCC. Climate Change 2014 impacts, adaptation, and vulnerability part A: global and sectoral aspects working group II contribution to the fifth assessment report of the intergovernmental panel on climate change foreword (Cambridge, New York, 2014).

18.
IPBES. The IPBES regional assessment report on biodiversity and ecosystem services for Europe and Central Asia (Secretariat of the Intergovernmental Science-Policy Platform on Biodiversity and Ecosystem Services, Bonn, 2018).

19.
Lomolino, M. V., Riddle, B. R. & Whittaker, R. J. Biogeography Biological Diversity Across Space and Time. 5th edn (Sinauer Associates, Sunderland, MA, 2017).

20.
Saladin, B. et al. Environment and evolutionary history shape phylogenetic turnover in European tetrapods. Nat. Commun. 10, 1–9 (2019).

21.
Eiserhardt, W. L., Svenning, J. C., Baker, W. J., Couvreur, T. L. P. & Balslev, H. Dispersal and niche evolution jointly shape the geographic turnover of phylogenetic clades across continents. Sci. Rep. 3, 1164 (2013).
CAS  PubMed  PubMed Central  Google Scholar 

22.
Rampal, E. et al. A minimal model for the latitudinal diversity gradient suggests a dominant role for ecological limits. Am. Nat. 194, E122–E133 (2019).

23.
Wiens, J. J. & Donoghue, M. J. Historical biogeography, ecology and species richness. Trends Ecol. Evol. 19, 639–644 (2004).
PubMed  Google Scholar 

24.
Baselga, A. Partitioning the turnover and nestedness components of beta diversity. Glob. Ecol. Biogeogr. 19, 134–143 (2010).
Google Scholar 

25.
Pointing, S. B., Bollard-Breen, B. & Gillman, L. N. Diverse cryptic refuges for life during glaciation. Proc. Natl Acad. Sci. USA 111, 5452–5453 (2014).
ADS  CAS  PubMed  Google Scholar 

26.
Keppel, G. et al. Refugia: identifying and understanding safe havens for biodiversity under climate change. Glob. Ecol. Biogeogr. 21, 393–404 (2012).
Google Scholar 

27.
Taberlet, P., Fumagalli, L., Wust-Saucy, A. G. & Cosson, J. F. Comparative phylogeography and postglacial colonization routes in Europe. Mol. Ecol. 7, 453–464 (1998).
CAS  PubMed  Google Scholar 

28.
Hewitt, G. M. Some genetic consequences of ice ages, and their role in divergence and speciation. Biol. J. Linn. Soc. 58, 247–276 (1996).
Google Scholar 

29.
Svenning, J. C. Deterministic Plio-Pleistocene extinctions in the European cool-temperate tree flora. Ecol. Lett. 6, 646–653 (2003).
Google Scholar 

30.
Huntley, B. & Birks, H. J. B. An Atlas of Past and Present Pollen Maps for Europe: 0–13000 Years Ago. (Cambridge University Press, Cambridge, 1983).

31.
Petit, R. J. et al. Glacial refugia: hotspots but not melting pots of genetic diversity. Science 300, 1563–1565 (2003).
ADS  CAS  PubMed  Google Scholar 

32.
Davis, M. B. & Shaw, R. G. Range shifts and adaptive responses to quaternary climate change. Science 292, 673–679 (2001).
ADS  CAS  PubMed  Google Scholar 

33.
Hewitt, G. M. Genetic consequences of climatic oscillations in the quaternary. Philos. T. R. Soc. B 359, 183–195 (2004).
CAS  Google Scholar 

34.
Tzedakis, P. C., Emerson, B. C. & Hewitt, G. M. Cryptic or mystic? Glacial tree refugia in Northern Europe. Trends Ecol. Evol. 28, 696–704 (2013).
CAS  PubMed  Google Scholar 

35.
Garcia Molinos, J., Schoeman, D. S., Brown, C. J. & Burrows, M. T. VoCC: an r package for calculating the velocity of climate change and related climatic metrics. Methods Ecol. Evol. 00, 1–8 (2019).
Google Scholar 

36.
Nobis, M. P. & Normand, S. KISSMig—a simple model for R to account for limited migration in analyses of species distributions. Ecography 37, 1282–1287 (2014).
Google Scholar 

37.
Thuiller, W., Lavorel, S., Araujo, M. B., Sykes, M. T. & Prentice, I. C. Climate change threats to plant diversity in Europe. Proc. Natl Acad. Sci. USA 102, 8245–8250 (2005).
ADS  CAS  PubMed  Google Scholar 

38.
Svenning, J. C. & Skov, F. Limited filling of the potential range in European tree species. Ecol. Lett. 7, 565–573 (2004).
Google Scholar 

39.
Svenning, J. C., Normand, S. & Skov, F. Range filling in European trees. J. Biogeogr. 33, 2018–2221 (2006).
Google Scholar 

40.
Svenning, J. C. & Skov, F. Could the tree diversity pattern in Europe be generated by postglacial dispersal limitation? Ecol. Lett. 10, 453–460 (2007).
PubMed  Google Scholar 

41.
Normand, S. et al. Postglacial migration supplements climate in determining plant species ranges in Europe. Proc. R. Soc. B 278, 3644–3653 (2011).
PubMed  Google Scholar 

42.
Arita, H. T. & Rodriguez, P. Geographic range, turnover rate and the scaling of species diversity. Ecography 25, 541–550 (2002).
Google Scholar 

43.
Slatyer, R. A., Hirst, M. & Sexton, J. P. Niche breadth predicts geographical range size: a general ecological pattern. Ecol. Lett. 16, 1104–1114 (2013).
PubMed  Google Scholar 

44.
Kambach, S. et al. Of niches and distributions: range size increases with niche breadth both globally and regionally but regional estimates poorly relate to global estimates. Ecography 42, 467–477 (2019).
Google Scholar 

45.
Warren, M. S. et al. Rapid responses of British butterflies to opposing forces of climate and habitat change. Nature 414, 65–69 (2001).
ADS  CAS  PubMed  Google Scholar 

46.
Hortal, J. et al. Ice age climate, evolutionary constraints and diversity patterns of European dung beetles. Ecol. Lett. 14, 741–748 (2011).
PubMed  Google Scholar 

47.
Pedreschi, D. et al. Challenging the European southern refugium hypothesis: species-specific structures versus general patterns of genetic diversity and differentiation among small mammals. Glob. Ecol. Biogeogr. 28, 262–274 (2019).
Google Scholar 

48.
Dullinger, S. et al. Post-glacial migration lag restricts range filling of plants in the European Alps. Glob. Ecol. Biogeogr. 21, 829–840 (2012).
Google Scholar 

49.
Kuhne, G., Kosuch, J., Hochkirch, A. & Schmitt, T. Extra-Mediterranean glacial refugia in a Mediterranean faunal element: the phylogeography of the chalk-hill blue Polyommatus coridon (Lepidoptera, Lycaenidae). Sci. Rep. 7, 43533 (2017).

50.
Svenning, J. C., Normand, S. & Kageyama, M. Glacial refugia of temperate trees in Europe: insights from species distribution modelling. J. Ecol. 96, 1117–1127 (2008).
Google Scholar 

51.
Willis, K. J., Rudner, E. & Sumegi, P. The full-glacial forests of central and southeastern Europe. Quat. Res. 53, 203–213 (2000).
Google Scholar 

52.
Stewart, J. R. & Lister, A. M. Cryptic northern refugia and the origins of the modern biota. Trends Ecol. Evol. 16, 608–613 (2001).
Google Scholar 

53.
Schonswetter, P., Stehlik, I., Holderegger, R. & Tribsch, A. Molecular evidence for glacial refugia of mountain plants in the European Alps. Mol. Ecol. 14, 3547–3555 (2005).
CAS  PubMed  Google Scholar 

54.
Abbott, R. J. et al. Molecular analysis of plant migration and refugia in the Arctic. Science 289, 1343–1346 (2000).
ADS  CAS  PubMed  Google Scholar 

55.
Stewart, J. R., Lister, A. M., Barnes, I. & Dalen, L. Refugia revisited: individualistic responses of species in space and time. Proc. R. Soc. B 277, 661–671 (2010).
PubMed  Google Scholar 

56.
Lenoir, J. et al. Cross-scale analysis of the region effect on vascular plant species diversity in Southern and Northern European mountain ranges. PLoS ONE 5, e15734 (2010).

57.
IPCC. Climate Change 2013: The physical science basis. Contribution of working group I to the fifth assessment report of the intergovernmental panel on climate change. (Cambridge University Press, Cambridge, New York, 2013).

58.
Jalas, J. & Suominen, J. Atlas Florae Europaeae, vol. 1–12. (The committee for mapping the flora of Europe & Societas Biologica Fennica Vanamo, Helsinki, 1999).

59.
Kurtto, A., Sennikov, A. N. & Lampinen, R. Atlas Florae Europaeae, vol. 13–16. (The committee for mapping the flora of Europe & Societas Biologica Fennica Vanamo, Helsinki, 2013).

60.
TPL. The Plant List Version 1.1 (The Plant List, 2013). http://www.theplantlist.org/.

61.
Cayuela, L., Stein, A. & Oksanen, J. Taxonstand: taxonomic standardization of plant species names. R package version 2.1. https://CRAN.R-project.org/package=Taxonstand (2017).

62.
R Core Team. R: A Language and Environment for Statistical Computing (R Foundation for Statistical Computing, Vienna, Austria, 2018). https://www.R-project.org/.

63.
Qian, H. & Jin, Y. An updated megaphylogeny of plants, a tool for generating plant phylogenies and an analysis of phylogenetic community structure. J. Plant Ecol. 9, 233–239 (2016).
Google Scholar 

64.
The information resource for euro-mediterranean plant diversity (Euro+Med PlantBase, 2006). http://ww2.bgbm.org/EuroPlusMed/.

65.
Smith, S. A. & Brown, J. W. Constructing a broadly inclusive seed plant phylogeny. Am. J. Bot. 105, 302–314 (2018).
PubMed  Google Scholar 

66.
Zanne, A. E. et al. Three keys to the radiation of angiosperms into freezing environments. Nature 506, 89–92 (2014).
ADS  CAS  PubMed  Google Scholar 

67.
Qian, H. & Zhang, J. Using an updated time-calibrated family-level phylogeny of seed plants to test for non-random patterns of life forms across the phylogeny. J. Syst. Evol. 52, 423–430 (2014).
Google Scholar 

68.
Webb, C. O., Ackerly, D. D. & Kembel, S. W. Phylocom: software for the analysis of phylogenetic community structure and trait evolution. Bioinformatics 24, 2098–2100 (2008).
CAS  PubMed  Google Scholar 

69.
Redelings, B. D. & Holder, M. T. A supertree pipeline for summarizing phylogenetic and taxonomic information for millions of species. Peerj 5, e3058 (2017).

70.
Synthesis release Open Tree of Life version 9.1 (Open Tree of Life, 2017). https://tree.opentreeoflife.org/about/synthesis-release/v9.1.

71.
Magallon, S., Gomez-Acevedo, S., Sanchez-Reyes, L. L. & Hernandez-Hernandez, T. A metacalibrated time-tree documents the early rise of flowering plant phylogenetic diversity. N. Phytol. 207, 437–453 (2015).
Google Scholar 

72.
Hijmans, R. J. Raster: geographic data analysis and modeling, v. 2.9-23. https://CRAN.R-project.org/package=raster (2019).

73.
Garcia Molinos, J., Schoeman, D. S., Brown, C. J. & Burrows, M. T. VoCC: the velocity of climate change and related climatic metrics. R package version 1.0.0. https://github.com/JorGarMol/VoCC (2019).

74.
Burrows, M. T. et al. The pace of shifting climate in marine and terrestrial ecosystems. Science 334, 652–655 (2011).
ADS  CAS  PubMed  Google Scholar 

75.
Maiorano, L. et al. Building the niche through time: using 13,000 years of data to predict the effects of climate change on three tree species in Europe. Glob. Ecol. Biogeogr. 22, 302–317 (2013).
Google Scholar 

76.
Owens, H. L. & Guralnick, R. P. ClimateStability: an R package to estimate climate stability from the time-slice climatologies. Biodivers. Inform. 14, 8–13 (2019).
Google Scholar 

77.
Zweig, M. H. & Campbell, G. Receiver-operating characteristic (Roc) plots—a fundamental evaluation tool in clinical medicine. Clin. Chem. 39, 561–577 (1993).
CAS  PubMed  Google Scholar 

78.
Leprieur, F. et al. Quantifying phylogenetic beta diversity: distinguishing between ‘true’ turnover of lineages and phylogenetic diversity gradients. PLoS ONE 7, e42760 (2012).
ADS  CAS  PubMed  PubMed Central  Google Scholar 

79.
Baselga, A. & Orme, C. D. L. Betapart: an R package for the study of beta diversity. Methods Ecol. Evol. 3, 808–812 (2012).
Google Scholar 

80.
Grömping, U. Relative importance for linear regression in R: the package relaimpo. J. Stat. Softw. 17, 1–27 (2006).
Google Scholar 

81.
EEA. Elevation map of Europe. https://www.eea.europa.eu/data-and-maps/data/digital-elevation-model-of-europe (2004).

82.
Paradis, E. & Schliep, K. ape 5.0: an environment for modern phylogenetics and evolutionary analyses in R. Bioinformatics 35, 526–528 (2019).
CAS  PubMed  Google Scholar 

83.
Wickham, H. stringr: simple, consistent wrappers for common string operations. R package version 1.3.1. https://CRAN.R-project.org/package=stringr (2018).

84.
Beygelzimer, A. et al. FNN: fast nearest neighbor search algorithms and applications. R package version 1.1.2.1. https://CRAN.R-project.org/package=FNN (2018).

85.
Bivand, R., Keitt, T. & Rowlingson, B. rgdal: bindings for the ‘Geospatial’ Data Abstraction Library. R package version 1.3-4. https://CRAN.R-project.org/package=rgdal (2018).

86.
Neuwirth, E. RColorBrewer: ColorBrewer palettes. R package version 1.1-2. https://CRAN.R-project.org/package=RColorBrewer (2014).

87.
Nychka, D., Furrer, R., Paige, J. & Sain, S. Fields: tools for spatial data. R package version 9.6. https://CRAN.R-project.org/package=fields (2017). More

Occurrence of the target OCs
The following individual OCs were found in fish samples with a detection frequency lower than 50%: i.e. p,p’-methoxychlor, heptachlor, heptachlor exo-epoxide, CB-66, 77, 81, 105, 114, 118, 123, 126, 128, 138, 153, 156, 157,167, 169, 170, and 180. These OCs are not discussed later. The concentrations of seven OC compounds in the muscle samples of 41 Nile tilapia and 34 Redbelly tilapia are shown in Supplementary Table S4. Since there were no significant interspecies differences between Nile tilapia and Redbelly tilapia (t-test, p = 0.16), the results of OCs analysis will be reported by tilapia genus in this study.
Median concentrations of OCPs and PCBs in tilapia samples from the main rivers system in the southern Guangxi are summarized in Table 1. PCBs and OCPs were detected in the muscle of all tilapia samples. DDTs were the predominant contaminant with a median concentration of 15.2 ng/g lw, and endosulfan was the second most common contaminant with a median concentration of 12.2 ng/g lw. PCBs, Drins, HCB, and HCHs concentrations in the fish were relatively low with median concentrations between 1.37 and 9.11 ng/g lw. The concentrations of the various OC compounds measured in the tilapia samples in this study were lower than those in tilapia collected from Guangdong province, China9, Africa10,11,12, Europe, and America13,14,15 (see Supplementary Table S5 online). This study showed that the main rivers in the southern Guangxi have low levels of OCs pollution, and the fish muscle contamination might be related to the low levels of pollution in the water and sediment. According to data from the National Bureau of statistics of China, the gross output value of industry and agriculture in Guangxi has been lower than that of other provinces or regions in China in the past few decades16. Therefore, the low levels of OCs pollution found in this study area are mainly the result of lower pollution input. In addition, most of the study area is located in the tropics, which have a relatively high perennial temperature. A warm climate is very conducive to enhance the metabolism rate of OCs by organisms17. The metabolism of organic pollutants by organisms occurs under the catalysis of a series of enzymes18,19. Factors affecting the enzymatic reaction, such as enzyme concentration and temperature, will affect the metabolism of OCs in organisms. Temperature also affects the air–water partitioning, which influences the volatilization of chemical pollutants from water20. Thus, dissolved chemical concentrations tend to be higher in cooler than in warmer waters21. In alignment with this supposition, Sobek et al. (2010) reported a largely reduced difference in bioaccumulation factor of PCBs between the Arctic and the temperate food webs, after adjustment for temperature effect22.
Table 1 Organochlorine concentration [median (range), ng/g lw] in the wild tilapia from the main rivers in Guangxi, South China.
Full size table

Distribution characteristics
Spatial distribution of OCs
The spatial distributions of seven OC compounds are presented in Fig. 1. The spatial distribution did not show a gradient in selected OCs concentrations. The spatial distribution of OCPs was under the double influence of a global distillation effect and the usage of OCPs23. Human activities can affect the distribution of OCPs in hilly areas24. However, there was no significant correlation between elevation and the residues of OCPs in this study (non-parametric test, p  > 0.05) (Fig. S2). Therefore, the distribution pattern of OCs in this study was hardly affected by global distillation. High levels of OCPs were found in TD and GG, where endosulfan or DDTs were the predominant contributors. Endosulfan is a cyclodiene pesticide extensively used throughout the world to control a wide variety of insects and mites23. Endosulfan levels were remarkably higher (10–411 times) in tilapia samples from TD than in samples from other sites. This observation was consistent with the fact that the local fruit and vegetable farming industry is the primary income source in the TD25. Therefore, we believe that the high levels of endosulfan in this study could be attributed to local pesticide practices specific to pest control needs over a short period26. Similarly, the higher levels of DDTs observed in GG also might be related to local short-term agricultural activities.
Figure 1

Spatial variations of log-transformed concentrations of OC compounds (ng/g lw) residues in wild tilapia from the main rivers in Guangxi, South China. TD: Tiandong County; LA: Longan County; CZ: Chongzuo City; FS: Fusui City; NN: Nanning City; GG: Guigang City; WX: Wuxuan County; PN: Pingnan City; TX: Tengxian County; WZ: Wuzhou City.

Full size image

PCBs are ubiquitous in tilapia samples from the study area, with a detection rate of 100%. In contrast to the OCP compounds, the overall trend of the PCBs was fairly homogenous. A relatively high median PCB concentration was detected in tilapia samples from TX, while slightly lower concentrations were detected from PN. There were no significant differences among different sampling locations (t-test, p  > 0.05). The minor differences could be explained by the migration and spread of PCBs in the environment. The limited historical use of PCBs in the present study area is another important factor contributing to this phenomenon25.
Spatial differences in pollutant metabolites
The ratio of parent compounds to their metabolites can provide useful information for the diagnosis of their sources23,24,27. The scatter plots for isomeric ratios of selected OCPs are shown in Fig. 2.
Figure 2

Scatter plots of molecular indices to identify DDTs (a) and endosulfan (b) contamination sources.

Full size image

The ratios between p,p’-DDT, p,p’-DDE and p,p’-DDD have been regarded as an indication of increasing or decreasing inputs to the environment. A ratio of (p,p’-DDE + p,p’-DDD)/p,p’-DDT  1.3) were found in two fish samples from FS (1.80) and CZ (1.71) districts, which indicates that Dicofol may be the main contributor to DDTs in these areas. In summary, the DDT residues in wild tilapia from rivers of the southern Guangxi originated mainly from the historical application of Dicofol and technical DDTs, whereas recent application of technical DDTs are indicated in TD.
Technical endosulfan includes two active ingredients: α-endosulfan (70%) and β-endosulfan (30%)23. Because α-endosulfan decomposes more easily than β-endosulfan, a α-/β-endosulfan ratio of  2.33) present in the tilapia samples from TD and FS, indicate continual use of endosulfan in these areas. In the other sites, those ratios are all below 2.33, indicating there was no recent application of technical endosulfan in that area. It is noteworthy that one sample (from TX site) contained β-endosulfan at a level below the limit of detection, but had appreciable levels of α-endosulfan, which may have been transported in from other areas. Because the Henry’s law constant for α-endosulfan is higher than the constant for β-endosulfan, there is a greater tendency for α-endosulfan to evaporate from the surface medium to air23,30.
The concentrations of ten PCB congeners in the present study area are illustrated in Fig. 3. Using degree of chlorination, these congeners can be divided into light PCBs (2–3 chlorines), medium PCBs (4–6 chlorines), and heavy PCBs (7–10 chlorines). The PCB sources of the 75 fish samples can be classified into the same categories since the PCBs in all sampling sites generally exhibited the following order: heavy PCBs (63.3–86.1% of ∑10PCBs)  > medium PCBs (9.72–18.2% of ∑10PCBs)  > light PCBs (4.66–18.3% of ∑10PCBs). The higher residual content of heavy PCBs may be related to historical production and use, or relate to their stably chemical structure31. Tri-CBs and penta-CBs were the major PCB products manufactured in China from 1965 until they were banned in 197431,32. The proportion of these compounds in PCBs was only 2.86–24.7% in this study. On the other hand, light PCBs have higher volatility and a lower octanol–water distribution coefficient than heavy PCBs33. Once absorbed into the organisms, light PCBs are usually more rapidly metabolized than the more highly chlorinated congeners34. Our results also indicated that the proportion of deca-CBs in heavy PCBs and PCBs was 78.3–98.1% and 51.8–88.6%, respectively. And the sampling sites with high deca-CBs ratio were distributed in the main agricultural farming areas (middle and upper reaches of rivers). And China banned the production of deca-CBs as early as 197435. Therefore, we believe that historical heritage was the main source of deca-CBs in the study area.
Figure 3

Composition profiles of PCB congeners in the main rivers from Guangxi, South China.

Full size image

Correlation among biological parameters and OC compounds
We studied the effects of biological parameters (including total length, body mass, age, and lipid content) on the bioaccumulation of individual contaminants in tissue samples of wild tilapia (based on dry weight). The loading and scores plot of PCA based on the concentrations of OCs in the tilapia muscle samples are displayed in Fig. 4. The PC1 explained 58.6% of the total variance and PC2 accounted for 23.5% of the variance. Table S7 lists the correlation coefficients between OC compounds and biological parameters, while the correlation coefficients between OC congeners and biological parameters are listed in Table S8. A significant relationship between growth parameters (i.e. total length, age, and body mass) was found in the tilapia samples, but only age and lipid content were significantly correlated (p  More

Environmental parameters of deep ocean seawater and sediment
Challenger Deep seawater and surface sediment samples were taken from its entire ~11,000 depth profile (Fig. 1a and Supplementary Table 1). The clines in temperature (29.8 °C in surface waters, decreasing to ~1.0 °C below 3000 m) and pressure (0.1 MPa in surface waters to ~104 MPa at the bottom of the trench) were recorded. The waters were oxic throughout the water column and the salinity ranged between 34 and 35 Practical Salinity Units (PSUs) (Supplementary Table 1). Seawater total DMSP and DMS concentrations were similar to those in previous studies21,26,27,28 and were positively correlated with Chl-a levels, being highest in the Chl-a maximum layer (10.51 × 10−3 nmol ml−1 DMSP and 4.97 × 10−3 nmol ml−1 DMS) and at lower but relatively stable levels (0.96–2.39 × 10−3 nmol ml−1 DMSP and 0.15–1.06 × 10−3 nmol ml−1 DMS) in the aphotic waters below 200 m (Fig. 1b, c and Supplementary Table 1). It should be noted that a small portion of this ‘background DMSP’ (3 μm) bacteria, which dominated the metagenomes of both these fractions, contained DMSP biosynthesis and catabolic genes (Figs. 2a and 3a, and Table 1), indicating that size fractionation is not a foolproof method of separating DMSP-producing bacteria from phytoplankton. Bacteria with dysB were shown by qPCR and metagenomics to be relatively abundant in the surface waters (dsyB total abundance of 2.61 × 105 copies L−1; 0.78–0.98% of surface water bacteria) representing ~3.49–4.38 × 103 bacteria ml−1 of surface seawater. These numbers are comparable to those predicted from the ocean microbial reference gene catalog metagenomic database (OM-RGC)37 in Williams et al.10, at ~4.8–9.6 × 103 bacteria ml−1. The abundance of these potential DMSP-producing bacteria initially decreased in 1000–2000 m deep seawater samples (3.46 × 104 copies L−1; ~0.43% bacteria at these depths), but then steadily increased with depth to reach maximal levels at 10,500 m (3.95 × 106 copies L−1; 4.03% of bacteria at 10,500 m), which were up to 15-fold higher than in the surface water (Fig. 2b, Table 1, and Supplementary Table 2). All detected dsyB sequences, including 37/162 metagenome assembled genomes (MAGs), were Alphaproteobacterial, mainly Rhodobacterales, Rhizobiales, and Rhodospirillales (Supplementary Data 1). At the genus level, Pseudooceanicola and Roseovarius were the most abundant potential DMSP producers at all depths, with much higher abundances (P  More

1.
Tilman, D. & Downing, J. A. Biodiversity and stability in grasslands. Nature 367, 363 (1994).
ADS  Google Scholar 
2.
Watkinson, A. R. & Ormerod, S. J. Grasslands, grazing and biodiversity: editors’ introduction. J. Appl. Ecol. 38, 233–237 (2001).
Google Scholar 

3.
Dengler, J., Janišová, M., Török, P. & Wellstein, C. Biodiversity of Palaearctic grasslands: a synthesis. Agric. Ecosyst. Environ. 182, 1–14 (2014).
Google Scholar 

4.
Dover, J. W., Spencer, S., Collins, S., Hadjigeorgiou, I. & Rescia, A. Grassland butterflies and low intensity farming in Europe. J. Insect Conserv. 15, 129–137 (2011).
Google Scholar 

5.
Morelli, F. High nature value farmland increases taxonomic diversity, functional richness and evolutionary uniqueness of bird communities. Ecol. Indic. 90, 540–546 (2018).
Google Scholar 

6.
Morelli, F., Benedetti, Y. & Tryjanowski, P. Introduction. In Birds as Useful Indicators of High Nature Value Farmlands (eds Morelli, F. & Tryjanowski, P.) 1–26 (Springer, Berlin, 2017). https://doi.org/10.1007/978-3-319-50284-7_1.
Google Scholar 

7.
Sutcliffe, L. M. E. et al. Harnessing the biodiversity value of Central and Eastern European farmland. Divers. Distrib. 21, 722–730 (2015).
Google Scholar 

8.
Donald, P. F., Green, R. E. & Heath, M. F. Agricultural intensification and the collapse of Europe’s farmland bird populations. Proc. R. Soc. Lond. Ser. B Biol. Sci. 268, 25–29 (2001).
Google Scholar 

9.
Donald, P. F., Pisano, G., Rayment, M. D. & Pain, D. J. The Common Agricultural Policy, EU enlargement and the conservation of Europe’s farmland birds. Agric. Ecosyst. Environ. 89, 167–182 (2002).
Google Scholar 

10.
Fragoso, R., Marques, C., Lucas, M. R., Martins, M. B. & Jorge, R. The economic effects of common agricultural policy on Mediterranean montado/dehesa ecosystem. J. Policy Model. 33, 311–327 (2011).
Google Scholar 

11.
Ribeiro, P. F. et al. Modelling farming system dynamics in high nature value farmland under policy change. Agric. Ecosyst. Environ. 183, 138–144 (2014).
Google Scholar 

12.
Suárez, F., Naveso, M. A. & De Juana, E. Farming in the drylands of Spain: birds of the pseudosteppes. In Farming and Birds in Europe. The common Agricultural Policy and its Implications for Bird Conservation (eds Pain, D. & Pienkowsky, M.) 297–330 (Academic Press, New York, 1997).
Google Scholar 

13.
Hoogeveen, Y., Petersen, J. E., Balazs, K. & Higuero, I. High Nature Value Farmland: Characteristics, Trends and Policy Challenges. EEA Report. No 1/2004. European Environment Agency, Copenhagen, Denmark (2004).

14.
Moreira, F., Pinto, M. J., Henriques, I. & Marques, T. The importance of low-intensive farming systems for fauna, flora and habitats protected under the european “birds” and “habitats” directives: is agriculture essential for preserving biodiversity in the mediterranean region? In Trends in Biodiversity Research (ed. Burk, A. R.) 117–145 (Nova Science Publishers, Huappauge, 2005).
Google Scholar 

15.
Lomba, A. et al. Mapping and monitoring high nature value farmlands: challenges in European landscapes. J. Environ. Manag. 143, 140–150 (2014).
Google Scholar 

16.
Delgado, A. & Moreira, F. Bird assemblages of an Iberian cereal steppe. Agric. Ecosyst. Environ. 78, 65–76 (2000).
Google Scholar 

17.
Ribeiro, P. F. et al. An applied farming systems approach to infer conservation-relevant agricultural practices for agri-environment policy design. Land Use Policy 58, 165–172 (2016).
Google Scholar 

18.
Stoate, C. et al. Ecological impacts of early 21st century agricultural change in Europe—a review. J. Environ. Manag. 91, 22–46 (2009).
CAS  Google Scholar 

19.
Faria, N., Morales, M. B. & Rabaça, J. E. Exploring nest destruction and bird mortality in mown Mediterranean dry grasslands: an increasing threat to grassland bird conservation. Eur. J. Wildl. Res. 62, 663–671 (2016).
Google Scholar 

20.
Santana, J. et al. Using beta diversity to inform agricultural policies and conservation actions on Mediterranean farmland. J. Appl. Ecol. 54, 1825–1835 (2017).
Google Scholar 

21.
Trombulak, S. C. & Frissell, C. A. Review of ecological effects of roads on terrestrial and aquatic communities. Conserv. Biol. 14, 18–30 (2000).
Google Scholar 

22.
Bernardino, J. et al. Bird collisions with power lines: state of the art and priority areas for research. Biol. Conserv. 222, 1–13 (2018).
Google Scholar 

23.
Loss, S. R., Will, T. & Marra, P. P. Direct mortality of birds from anthropogenic causes. Annu. Rev. Ecol. Evol. Syst. 46, 99–120 (2015).
Google Scholar 

24.
Hernández-Matías, A., Real, J., Parés, F. & Pradel, R. Electrocution threatens the viability of populations of the endangered Bonelli’s eagle (Aquila fasciata) in Southern Europe. Biol. Conserv. 191, 110–116 (2015).
Google Scholar 

25.
Shaw, J. M., Reid, T. A., Schutgens, M., Jenkins, A. R. & Ryan, P. G. High power line collision mortality of threatened bustards at a regional scale in the Karoo, South Africa. Ibis (Lond. 1859) https://doi.org/10.1111/ibi.12553 (2017).
Article  Google Scholar 

26.
Borda-de-Água, L., Grilo, C. & Pereira, H. M. Modeling the impact of road mortality on barn owl (Tyto alba) populations using age-structured models. Ecol. Model. 276, 29–37 (2014).
Google Scholar 

27.
Reijnen, R., Foppen, R. & Meeuwsen, H. The effects of traffic on the density of breeding birds in Dutch agricultural grasslands. Biol. Conserv. 75, 255–260 (1996).
Google Scholar 

28.
Mcnew, L. B., Hunt, L. M., Gregory, A. J., Wisely, S. M. & Sandercock, B. K. Effects of wind energy development on nesting ecology of greater prairie-chickens in fragmented grasslands. Conserv. Biol. 28, 1089–1099 (2014).
PubMed  PubMed Central  Google Scholar 

29.
Wolfe, D. H., Patten, M. A., Shochat, E., Pruett, C. L. & Sherrod, S. K. Causes and patterns of mortality in lesser prairie-chickens Tympanuchus pallidicinctus and implications for management. Wildl. Biol. 13, 95–104 (2007).
Google Scholar 

30.
Shaw, J. M., Jenkins, A. R., Smallie, J. J. & Ryan, P. G. Modelling power-line collision risk for the Blue Crane Anthropoides paradiseus in South Africa. Ibis (Lond. 1859) 152, 590–599 (2010).
Google Scholar 

31.
Birdlife International. The IUCN Red List of Threatened Species 2018 (2018). www.iucnredlist.org. Accessed 2nd August 2019.

32.
Faria, N. Implications of Dry Grassland Management in the Ecology and Conservation of Grassland Birds in South Portugal (Universidad Autónoma de Madrid, Madrid, 2015).
Google Scholar 

33.
Iñigo, A. & Barov, B. Action plan for the Little Bustard Tetrax tetrax in the European Union. Report. SEO| BirdLife and BirdLife International for the European Commission (2010).

34.
Morales, M. B., García, J. T. & Arroyo, B. Can landscape composition changes predict spatial and annual variation of little bustard male abundance?. Anim. Conserv. 8, 167–174 (2005).
Google Scholar 

35.
Moreira, F. et al. Mosaic-level inference of the impact of land cover changes in agricultural landscapes on biodiversity: a case-study with a threatened grassland bird. PLoS ONE 7, e38876 (2012).
ADS  MathSciNet  CAS  PubMed  PubMed Central  Google Scholar 

36.
Silva, J. P., Palmeirim, J. M. & Moreira, F. Higher breeding densities of the threatened little bustard Tetrax tetrax occur in larger grassland fields: implications for conservation. Biol. Conserv. 143, 2553–2558 (2010).
Google Scholar 

37.
Silva, J. P. et al. EU protected area network did not prevent a country wide population decline in a threatened grassland bird. PeerJ 6, e4284 (2018).
PubMed  PubMed Central  Google Scholar 

38.
García de la Morena, Bota, G., Mañosa, S. & Morales, M. B. El sisón común en España. II Censo Nacional (2016). Report (2018).

39.
Traba, J. & Morales, M. B. The decline of farmland birds in Spain is strongly associated to the loss of fallowland. Sci. Rep. 9, 9473 (2019).
ADS  PubMed  PubMed Central  Google Scholar 

40.
Marcelino, J. et al. Tracking data of the Little Bustard Tetrax tetrax in Iberia shows high anthropogenic mortality. Bird Conserv. Int. https://doi.org/10.1017/S095927091700051X (2017).
Article  Google Scholar 

41.
Bevanger, K. Biological and conservation aspects of bird mortality caused by electricity power lines: a review. Biol. Conserv. 86, 67–76 (1998).
Google Scholar 

42.
Janss, G. F. E. Avian mortality from power lines: a morphologic approach of a species-specific mortality. Biol. Conserv. 95, 353–359 (2000).
Google Scholar 

43.
Martin, G. R. Understanding bird collisions with man-made objects: a sensory ecology approach. Ibis. 153, 239–254 (2011).
Google Scholar 

44.
Martin, G. R. & Shaw, J. M. Bird collisions with power lines: failing to see the way ahead?. Biol. Conserv. 143, 2695–2702 (2010).
Google Scholar 

45.
Marques, A. T., Martins, R. C., Silva, J. P., Palmeirim, J. M. & Moreira, F. Power line routing and configuration as major drivers of collision risk in two bustard species. Oryx https://doi.org/10.1017/S0030605319000292 (2020).
Article  Google Scholar 

46.
Silva, J. P. et al. A spatially explicit approach to assess the collision risk between birds and overhead power lines: a case study with the little bustard. Biol. Conserv. 170, 256–263 (2014).
Google Scholar 

47.
García, J., Suárez-Seoane, S., Miguélez, D., Osborne, P. E. & Zumalacárregui, C. Spatial analysis of habitat quality in a fragmented population of little bustard (Tetrax tetrax): implications for conservation. Biol. Conserv. 137, 45–56 (2007).
Google Scholar 

48.
Osborne, P. E. & Suárez-Seoane, S. Identifying core areas in a species’ range using temporal suitability analysis: an example using little bustards Tetrax Tetrax L. in Spain. Biodivers. Conserv. 16, 3505–3518 (2007).
Google Scholar 

49.
Santangeli, A. & Dolman, P. M. Density and habitat preferences of male little bustard across contrasting agro-pastoral landscapes in Sardinia (Italy). Eur. J. Wildl. Res. 57, 805–815 (2011).
Google Scholar 

50.
Santos, M. et al. Impacts of land use and infrastructural changes on threatened Little Bustard Tetrax tetrax breeding populations: quantitative assessments using a recently developed spatially explicit dynamic modelling framework. Bird Conserv. Int. 26, 418–435 (2016).
Google Scholar 

51.
Suárez-Seoane, S., Osborne, P. E. & Alonso, J. C. Large-scale habitat selection by agricultural steppe birds in Spain: identifying species-habitat responses using generalized additive models. J. Appl. Ecol. 39, 755–771 (2002).
Google Scholar 

52.
Silva, J. P. et al. Estimating the influence of overhead transmission power lines and landscape context on the density of little bustard Tetrax tetrax breeding populations. Ecol. Model. 221, 1954–1963 (2010).
Google Scholar 

53.
Morales, M. B., Traba, J., Carriles, E., Delgado, M. P. & de la Morena, E. L. G. Sexual differences in microhabitat selection of breeding little bustards Tetrax tetrax: ecological segregation based on vegetation structure. Acta Oecologica 34, 345–353 (2008).
ADS  Google Scholar 

54.
Faria, N., Rabaça, J. E. & Morales, M. B. The importance of grazing regime in the provision of breeding habitat for grassland birds: the case of the endangered little bustard (Tetrax tetrax). J. Nat. Conserv. 20, 211–218 (2012).
Google Scholar 

55.
Silva, J. P., Estanque, B., Moreira, F. & Palmeirim, J. M. Population density and use of grasslands by female Little Bustards during lek attendance, nesting and brood-rearing. J. Ornithol. 155, 53–63 (2014).
Google Scholar 

56.
INE. Statistical data: Database. (2019). https://www.ine.pt/xportal/xmain?xpid=INE&xpgid=ine_base_dados. Accessed 9th May 2019.

57.
Gameiro, J., Silva, J. P., Franco, A. M. A. & Palmeirim, J. M. Effectiveness of the European Natura 2000 network at protecting Western Europe’s agro-steppes. Biol. Conserv. 248, 108681 (2020).
Google Scholar 

58.
Beja, P. et al. Predators and livestock reduce bird nest survival in intensive Mediterranean farmland. Eur. J. Wildl. Res. 60, 249–258 (2014).
Google Scholar 

59.
van der Wal, R. & Palmer, S. C. Is breeding of farmland wading birds depressed by a combination of predator abundance and grazing?. Biol. Lett. 4, 256–258 (2008).
PubMed  PubMed Central  Google Scholar 

60.
Lane, S. J., Alonso, J. C. & Martín, C. A. Habitat preferences of great bustard Otis tarda flocks in the arable steppes of central Spain: are potentially suitable areas unoccupied?. J. Appl. Ecol. 38, 193–203 (2001).
Google Scholar 

61.
Ahlering, M. A., Johnson, D. H. & Faaborg, J. Conspecific attraction in a grassland bird, the Baird’s Sparrow. J. Field Ornithol. 77, 365–371 (2006).
Google Scholar 

62.
Tarjuelo, R. et al. Not only habitat but also sex: factors affecting spatial distribution of Little Bustard Tetrax tetrax families. Acta Ornithol. 48, 119–128 (2013).
Google Scholar 

63.
Reino, L. et al. Effects of changed grazing regimes and habitat fragmentation on Mediterranean grassland birds. Agric. Ecosyst. Environ. 138, 27–34 (2010).
Google Scholar 

64.
Walters, K., Kosciuch, K. & Jones, J. Can the effect of tall structures on birds be isolated from other aspects of development?. Wildl. Soc. Bull. 38, 250–256 (2014).
Google Scholar 

65.
Fahrig, L. & Rytwinski, T. Effects of roads on animal abundance: an empirical review and synthesis. Ecol. Soc. 14, 21 (2009).
Google Scholar 

66.
Tryjanowski, P. et al. Conservation of farmland birds faces different challenges in Western and Central-Eastern Europe. Acta Ornithol. 46, 1–12 (2011).
Google Scholar 

67.
Gudka, M., Santos, C. D., Dolman, P. M., Abad-Gómez, J. M. & Silva, J. P. Feeling the heat: elevated temperature affects male display activity of a lekking grassland bird. PLoS ONE 14, e0221999 (2019).
CAS  PubMed  PubMed Central  Google Scholar 

68.
Silva, J. P., Catry, I., Palmeirim, J. M. & Moreira, F. Freezing heat: thermally imposed constraints on the daily activity patterns of a free-ranging grassland bird. Ecosphere 6, art119 (2015).

69.
Alonso, H. et al. Male post-breeding movements and stopover habitat selection of an endangered short-distance migrant, the Little Bustard Tetrax tetrax. Ibis (Lond. 1859) 162, 279–292 (2020).
Google Scholar 

70.
García de la Morena, E. L. et al. Migration patterns of Iberian little bustards Tetrax tetrax. Ardeola 62, 95–112 (2015).
Google Scholar 

71.
Silva, J. P., Faria, N. & Catry, T. Summer habitat selection and abundance of the threatened little bustard in Iberian agricultural landscapes. Biol. Conserv. 139, 186–194 (2007).
Google Scholar 

72.
Equipa Atlas. Atlas das aves nidificantes em Portugal (1999–2005). (Instituto da Conservaçăo da Natureza e da Biodiversidade, Sociedade Portuguesa para o Estudo das Aves, Parque Natural da Madeira e Secretaria Regional do Ambiente e do Mar, 2008).

73.
De Juana, E. & Martínez, C. Distribution and conservation status of Little bustard Tetrax tetrax in the Iberian Peninsula. Ardeola 43, 157–167 (1996).
Google Scholar 

74.
Delgado, A. & Moreira, F. Between-year variations in Little Bustard Tetrax tetrax population densities are influenced by agricultural intensification and rainfall. Ibis (Lond. 1859) 152, 633–642 (2010).
Google Scholar 

75.
DGT. Especificações técnicas da Carta de Uso e Ocupação do Solo de Portugal Continental para 1995, 2007, 2010 e 2015. 103 (2018).

76.
INE. Recenseamento Agrícola 1999—Análise de resultados (2001).

77.
INE. Recenseamento Agrícola 2009—Análise dos principais resultados (2011).

78.
Haklay, M. & Weber, P. Openstreetmap: user-generated street maps. IEEE Pervasive Comput. 7, 12–18 (2008).
Google Scholar 

79.
R Core Team. R: a language and environment for statistical computing (2016).

80.
Bates, D., Mächler, M., Bolker, B. & Walker, S. Fitting linear mixed-effects models using lme4. J. Stat. Softw. 67, 1–48 (2015).
Google Scholar 

81.
Kuznetsova, A., Brockhoff, P. B. & Christensen, R. H. B. lmerTest package: tests in linear mixed effects models. J. Stat. Softw. 82, 1–26 (2017).
Google Scholar 

82.
Wood, S. N. Generalized Additive Models: An Introduction with R 2nd edn. (Chapman and Hall/CRC, Boca Raton, 2017).
Google Scholar 

83.
Zuur, A. F., Ieno, E. N., Walker, N., Saveliev, A. A. & Smith, G. M. Mixed Effects Models and Extensions in Ecology with R (Springer, New York, 2009). https://doi.org/10.1007/978-0-387-87458-6.
Google Scholar 

84.
Wood, S. N. mgcv: mixed GAM computation vehicle with automatic smoothness estimation. R package version 1.8-24. 302 (2018).

85.
Bjørnstad, O. N. & Falck, W. Nonparametric spatial covariance functions: estimation and testing. Environ. Ecol. Stat. 8, 53–70 (2001).
MathSciNet  Google Scholar 

86.
Rhodes, J. R., McAlpine, C. A., Zuur, A. F., Smith, G. M. & Ieno, E. N. Mixed Effects Models and Extensions in Ecology with R 469–492 (Springer, Berlin, 2009). https://doi.org/10.1007/978-0-387-87458-6_21.
Google Scholar 

87.
Bjørnstad, O. N. ncf: spatial nonparametric covariance functions. R package version 1.1–7. (2016).

88.
QGIS Development Team. QGIS Geographic Information System. Open Source Geospatial Foundation Project. https://qgis.osgeo.org (2017). More

1.
Schuster, P. F. et al. Permafrost stores a globally significant amount of mercury. Geophys. Res. Lett. 45, 1463–1471 (2018).
ADS  CAS  Article  Google Scholar 
2.
Smith-Downey, N. V., Sunderland, E. M. & Jacob, D. J. Anthropogenic impacts on global storage and emissions of mercury from terrestrial soils: Insights from a new global model. J. Geophys. Res. 115, G03008 (2010).
ADS  Article  Google Scholar 

3.
Zimov, S. A. et al. Permafrost carbon: Stock and decomposability of a globally significant carbon pool. Geophys. Res. Lett. 33, https://doi.org/10.1029/2006GL027484 (2006).

4.
Romanovsky, V., Grosse, G. & Marchenko, S. Past, present and future of permafrost in a changing world. Geo. Soc. Am. 40, 397 (2008).
Google Scholar 

5.
Biskaborn et al. Permafrost is warming at a global scale. Nat. Comm. 10, 264 (2019).
ADS  Article  Google Scholar 

6.
Koven, C. D., Riley, W. J. & Stern, A. Analysis of permafrost thermal dynamics and response to climate change in the CMIP5 Earth System Models. J. Clim. V26, 1887–1900 (2013).
ADS  Google Scholar 

7.
McGuire, A. D. et al. Dependence of the evolution of carbon dynamics in the northern permafrost region on the trajectory of climate change. Proc. Natl Acad. Sci. USA 115/15, 3882–3887 (2018).
ADS  Article  Google Scholar 

8.
Driscoll, C. T., Mason, R. P., Chan, H. M., Jacob, D. J. & Pirrone, N. Mercury as a global pollutant: sources, pathways, and effects. Environ. Sci. Technol. 47, 4967–4983 (2013).
ADS  CAS  Article  Google Scholar 

9.
Obrist, D. et al. Tundra uptake of atmospheric elemental mercury drives Arctic mercury pollution. Nature 547, 201–204 (2017).
ADS  CAS  Article  Google Scholar 

10.
Skyllberg, U., Bloom, P. R., Qian, J., Lin, C. M. & Bleam, W. F. Complexation of mercury(II) in soil organic matter: EXAFS evidence for linear two-coordination with reduced sulfur groups. Environ. Sci. Technol. 40, 4174–4180 (2006).
ADS  CAS  Article  Google Scholar 

11.
Giesler, R., Clemmensen, K. E., Wardle, D. A., Klaminder, J. & Bindler, R. Boreal forests sequester large amounts of mercury over millennial time scales in the absence of wildfire. Environ. Sci. Technol. 51, 2621–2627 (2017).
ADS  CAS  Article  Google Scholar 

12.
Arnold, J., Gustin, M. S. & Weisberg, P. J. Evidence for nonstomatal uptake of Hg by Aspen and translocation of Hg from foliage to tree rings in Austrian pine. Environ. Sci. Technol. 52, 1174–1182 (2018).
ADS  CAS  Article  Google Scholar 

13.
Clackett, S. P., Porter, T. J. & Lehnherr, I. 400-year record of atmospheric mercury from tree-rings in Northwestern Canada. Environ. Sci. Technol. 52, 9625–9633 (2018).
ADS  CAS  Article  Google Scholar 

14.
Lindberg, S. E., Hanson, P. J., Meyers, T. P. & Kim, K. H. Air/surface exchange of mercury vapor over forests – the need for a reassessment of continental biogenic emissions. Atm. Environ. 32, 895–908 (1998).
CAS  Article  Google Scholar 

15.
Jiskra, M. et al. Mercury deposition and re-emission pathways in boreal forest soils investigated with Hg isotope signatures. Environ. Sci. Technol. 49, 7188–7196 (2015).
ADS  CAS  Article  Google Scholar 

16.
Schuster, P. F. et al. Mercury export from the Yukon River Basin and potential response to a changing climate. Environ. Sci. Technol. 45, 9262–9267 (2011).
ADS  CAS  Article  Google Scholar 

17.
Schaefer, K. et al. Combined Simple Biosphere/Carnegie-Ames-Stanford Approach terrestrial carbon cycle model. J. Geophys. Res. 113, G03034 (2008).
Article  Google Scholar 

18.
Olson, C., Jiskra, M., Biester, H., Chow, J. & Obrist, D. Mercury in active-layer Tundra soils of Alaska: concentrations, pools, origins, and spatial distribution. Glob. Biogeochemical Cycles 32, 1058–1073 (2018).
ADS  CAS  Article  Google Scholar 

19.
Mikan, C. J., Schimel, J. P. & Doyle, A. P. Temperature controls of microbial respiration in arctic tundra soils above and below freezing. Soil Biol. Biochem. 34, 1785–1795 (2002).
CAS  Article  Google Scholar 

20.
Wang, Z. & Roulet, N. Comparison of plant litter and peat decomposition changes with permafrost thaw in a subarctic peatland. Plant Soil 417, 197–216 (2017).
CAS  Article  Google Scholar 

21.
Wickland, K. P. et al. Dissolved organic carbon and nitrogen release from boreal Holocene permafrost and seasonally frozen soils of Alaska. Environ. Res. Lett. 13, 065011 (2018).
ADS  Article  Google Scholar 

22.
Striegl, R. G., Aiken, G. R., Dornblaser, M. M., Raymond, P. A. & Wickland, K. P. A decrease in discharge-normalized DOC export by the Yukon River during summer through autumn. Geophys. Res. Lett. 32, L21413 (2005).
ADS  Article  Google Scholar 

23.
Walvoord, M. A. & Striegl, R. G. Increased groundwater to stream discharge from permafrost thawing in the Yukon River basin: potential impacts on lateral export of carbon and nitrogen. Geophys. Res. Lett. 34, L12402 (2007).
ADS  Article  Google Scholar 

24.
Holmes, C. D. et al. Global atmospheric model for mercury including oxidation by bromine atoms. Atmos. Chem. Phys. 10, 12037–12057 (2010).
ADS  CAS  Article  Google Scholar 

25.
Pacyna, J. M. et al. Current and future levels of mercury atmospheric pollution on a global scale. Atmos. Chem. Phys. 16, 12495–12511 (2016).
ADS  CAS  Article  Google Scholar 

26.
Schaefer, K., Zhang, T., Bruhwiler, L. & Barrett, A. P. Amount and timing of permafrost carbon release in response to climate warming. Tellus Series B Chem. Phys. Met. https://doi.org/10.1111/j1600-0889201100527x (2011).

27.
St Pierre, K. A. et al. Unprecedented increases in total and methyl mercury concentrations downstream of retrogressive thaw slumps in the western Canadian arctic. Environ. Sci. Technol. 52, 14099–14109 (2018).
ADS  Article  Google Scholar 

28.
EPA. Ambient water quality criteria for mercury. U.S. Environmental Protection Agency, 440/5-84-026 (https://www.epa.gov/sites/production/files/2019-03/documents/ambient-wqc-mercury-1984.pdf) (1984).

29.
Brumbaugh, W. G., Krabbenhoft, D. P., Helsel, D. R., Wiener, J. G., & Echols, K. R. A national pilot study of mercury contamination of aquatic ecosystems along multiple gradients: bioaccumulation in fish, Biological Science Report. USGS/BRD/BSR-2001-0009 (2001).

30.
Scudder, E. et al. Optimizing fish sampling for fish-mercury bioaccumulation factors. Chemosphere 135, 467–473 (2015).
ADS  Article  Google Scholar 

31.
National Research Council. Toxicological effects of methylmercury. https://doi.org/10.17226/9899 (The National Academies Press, Washington, DC, 2000).

32.
Borum, D., Manibusan, M. K., Schoeny, R., Winchester, E. L. Water quality criterion for the protection of human health: methylmercury. EPA-823-R-01-001 (U.S. Environmental Protection Agency, Washington, DC 20460, 2001).

33.
Schaefer, K. et al. Improving simulated soil temperatures and soil freeze/thaw at high-latitude regions in the Simple Biosphere/Carnegie-Ames-Stanford Approach model. J. Geophys. Res. 114, F02021 (2009).
ADS  Article  Google Scholar 

34.
Schaefer, K., Lantuit, H., Romanovsky, V. E., Schuur, E. A. G. & Witt, R. The impact of the permafrost carbon feedback on global climate. Env. Res. Lett. 9, 085003 (2014).
Article  Google Scholar 

35.
Schaefer, K. & Jafarov, E. A parameterization of respiration in frozen soils based on substrate availability. Biogeosciences 13, 1991–2001. www.biogeosciences.net/13/1991/2016/ (2016).

36.
Jafarov, E. & Schaefer, K. The importance of a surface organic layer in simulating permafrost thermal and carbon dynamics. Cryosphere 10, 465–475 (2016).
ADS  Article  Google Scholar 

37.
USGS, United States Geological Survey. Data inventory page for site 15565447-Yukon River at Pilot Station, Alaska. U.S. Geological Survey, https://waterdata.usgs.gov/nwis/inventory/site_no=15565447 (2019).

38.
Agnan, Y., Le Dantec, T., Moore, C. W., Edwards, G. C. & Obrist, D. New constraints on terrestrial surface atmosphere fluxes of gaseous elemental mercury using a global database. Environ. Sci. Technol. 50, 507–524 (2016).
ADS  CAS  Article  Google Scholar  More

Locating and collecting badgers
Members of the public, local authorities and countryside organisations in Wales reported the locations of found dead badgers to the Welsh Government or APHA Field Services who recorded the map reference details. The collection by APHA staff of badger carcases that met pre-defined criteria took place between 1st September 2014 and 31st December 2016. There were instances where the collection of reported carcases was not attempted. Reasons for not attempting collection included safety concerns arising from the specific location of the carcass or the non-availability of staff (or other necessary resource) to undertake the collection (supplementary table online ST1). In some further instances collection was attempted but was unsuccessful because either the reported carcase could not be found or the condition of the carcase made it unsuitable for further investigation (for example viscera were herniated externally through wounds, there was severe myiasis (flystrike), the carcase was distended with gas or it was flattened).
Depending on where they were found, the carcases were delivered to the APHA Veterinary Investigation Centres in Carmarthen and Shrewsbury and to the Wales Veterinary Science Centre in Aberystwyth (from 28/01/2016) where they were stored at between 2 and 8 °C for no more than four days before post-mortem examination.
Post-mortem examination and sampling
Of the 1863 dead badgers reported, 841 were collected and 681 (37% of reported carcasses) were suitable for post-mortem examination. The prevalence of bTB was calculated for suitable carcasses only. The sampling protocol was adapted from the standard and detailed protocols described in a comparison by Crawshaw and others24, so that fewer overall samples were taken than the detailed protocol but that the samples chosen were the ones most likely to yield M. bovis if present.
The initial external examination comprised the following: measuring the length from nose tip to tail base (cm), assessing and recording tooth wear, recording the sex of the animal and whether female animals were lactating, andrecording any evidence of vaccination: To temporarily identify vaccinated badgers, the guard hairs of the dorsal back (usually caudal) were trimmed and a coloured marker applied at the same site at the time of vaccination. During the external examination of the badger carcases, to attempt to detect recent vaccination, the skin and hair of the back was visually examined for guard hair trimming and coloured marker. Furthermore the examination recorded any evidence of trap injury, or of illegal interference with the animal and the presence and location of bite wounds.
A detailed examination of the lungs, pericardial sac, liver and kidneys was conducted at post-mortem examination. The lungs were examined by making multiple longitudinal incisions approximately one centimetre apart. At least four slices were made in the liver and three slices in the kidney24.
Each lymph node of all suitable badgers was incised at least once and a pool of lymph nodes (pool 1) consisting of the retropharyngeal, bronchial, mediastinal and hepatic lymph nodes (or as many as were detectable) was collected for subsequent bacteriological culture. If any gross internal lesions suggestive of tuberculosis were observed or if bite wounds were detected, the lesioned tissue and/or excised bite wounds were added to a separate container (pool 2). The pooled samples for bacteriological culture were preserved in 15 ml of 1% aqueous cetylpyridinium chloride. Samples were sent to the Animal Plant & Health Agency Laboratory inStarcross, Devon on the day of examination, for next day receipt. After taking samples, two or three incisions were made in the muscles of the anterior thigh of both hind legs to look for any potential adverse reactions to BCG vaccination24.
Culture and molecular typing
On receipt at APHA Starcross, the samples were washed in sterile 0.85% saline solution, homogenized by standard methods, inoculated onto six modified Middlebrook 7H11 agar slopes26 and incubated at 37 °C. Pool 1, and Pool 2 (if collected), were cultured separately. The slopes were examined weekly from the end of week 2 for a maximum of 12 weeks. M. bovis isolates were harvested when growth was sufficient for genotyping and sent to APHA Weybridge.
Genotyping was performed using spoligotyping27 and VNTR typing (Exact Tandem Repeat Loci A to F)9,28. Spoligotyping confirmed that the isolates were M. bovis. Genotypes of M. bovis were labelled according to the current APHA convention, using numbers to represent spoligotypes and lower case letters to represent the VNTR pattern within each spoligotype. The same genotyping methods were applied to the cultures of M. bovis from cattle slaughtered as part of the national bTB control programme. The cattle genotype home ranges were determined for 2015 (Data Systems Group, Dept of Epidemiological Sciences, Animal & Plant Health Agency Weybridge); for inclusion a genotype had to have been present for at least three years, on at least two different 5 km × 5 km grid squares, in the last 5 years (with a 10 km buffer applied).
Disease status in cattle herds
The disease status in cattle herds from the five TB regions of Wales has been calculated for all Spatial Units (for definition please refer to introduction) with at least 10 badger submissions and compared with prevalence estimates for the badger populations in these areas. The metric used to describe disease status in cattle is average annual herd prevalence during most of the AWBFD sampling period (2015–2016, number of herds under restriction at any point of year/registered live herds in region). Herd size is a known risk factor for bovine TB29; therefore, direct standardisation was used to adjust for varying herd sizes in the Spatial Units3,30,31. Briefly, annual stratum-specific prevalence was calculated for each Spatial Unit across four strata (reduced from six in the cited studies due to the small herd denominator in Wales Spatial Units) of herd size and then applied to the reference population, which comprised the sum of cattle populations across all Spatial Units. The standardised population was then used for herd-level disease measures, resulting in a standardised herd prevalence (Fig. 4).
Badger prevalence mapping
For the Spatial Unit badger prevalence map each Spatial Unit was labelled with the number of submissions, the number of positive results and the resulting prevalence estimate. When added to ArcMap (Esri ArcGIS 10.2.2) the Spatial Unit layer was symbolised using the prevalence value and six pre-defined range values or classes were applied and colour ramped. Spatial Units with less than 10 AWBFD submissions (“insufficient data”) were not colour ramped. All Spatial Units were labelled with the number of positives/number of submissions.
Statistical analysis
The prevalence of bTB in badgers was estimated among the sampled badgers as in previous studies3,8 with the underlying assumption that the carcasses collected were representative of the overall population.
Analysis was performed to test the null hypotheses (Ho) that:

There are no differences in badger bTB prevalence between the five TB Areas in Wales.

There was no change in overall badger bTB prevalence in Wales between the surveys in 2005–2006 and in 2014–2016.

There was no change in badger bTB prevalence within the five TB Areas between the surveys in 2005–2006 and in 2014–2016.

There was no change in badger bTB prevalence between the surveys in 2005–2006 and in 2014–2016 within the Intensive Action Area (IAA), site of a 2012–2015 badger vaccination trial.

There is no correlation between cattle herd prevalence and estimated bTB prevalence in dead badgers in different geographic regions.

There is no difference in cattle herd prevalence between areas with infected badgers and those with no evidence of bTB in badgers.

For statistical purposes, for the comparisons of prevalence estimates between the TB Areas and over time, both Intermediate TB Areas were combined. All data were tested for normality with the Kolmogorov–Smirnov test (SPSS Version 21 for Windows). Since the dependent variable in these analyses, prevalence, is a rate and fulfilled the criteria for normality, the z-test for comparisons of population proportions was used to test for the statistical significance of differences between prevalence estimates between the TB Areas and over time. A condition of the z-test is that each sample contains at least 10 observations in each category of the dependent variable and for comparisons between samples with less than 10 submissions in at least one category, Fishers Exact test was used32. To explore the correlation between prevalence in badgers and cattle, linear regression analysis was undertaken to calculate Pearson’s coefficient of correlation (SPPS 22 for Windows). Spatial Units with  > 10 badger submissions with one or more positives were compared with those which had none, using the two sample t test with unequal variances. In order to compare genotypes between badgers and cattle, as in the previous Wales badger survey3, the authors calculated the associations between the frequency distribution of the genotypes in badgers and the frequency distribution in cattle for each TB Area. In order to prevent ties and to account for the small number of positive badger submissions, frequencies were adjusted, by replacing them with their deviations from expected values which were calculated as (TB Area subtotal) × (Wales genotype subtotal)/(Wales grand total). A Spearman rank correlation between the ranks of genotypes in badgers and their ranks in cattle was then calculated for each TB Area. More