Ecology

Rhizosphere soil samplingA total of 20 rhizosphere soil samples (20 tomato plants) were collected at the flowering stage from a tomato field located in Qilin town, Jiangsu province, China, 118°57’ E, 32°03’ N, which had been infested by the pathogen Ralstonia solanacearum for more than 15 years [8]. After uprooting plants, excess soil was first gently shaken from the roots, and the remaining soil attached to roots was considered as rhizosphere soil. Each rhizosphere soil sample was then used for bacterial strain isolation.Isolation and identification of rhizobacteriaIsolationA total of 640 bacterial strains were isolated from the fresh rhizosphere soil samples, according to a previously established protocol [11]. Briefly, 1 g of each rhizosphere sample was mixed with 9 mL MS buffer solution (50 mM Tris-HCl [pH 7.5], 100 mM NaCl, 10 mM MgSO4, 0.01% gelatin) in a rotary shaker at 170 rpm min−1 for 30 min at 30 °C. After serial dilution in MS buffer solution, 100-μl volumes of the diluted soil suspensions were plated on 1/10 tryptone soy agar (1/10 TSA, 1.5 g L−1 tryptone, 0.5 g L−1 soytone, 0.5 g L−1 sodium chloride, and 15 g L−1 agar, pH 7.0). After a 48-h incubation at 30 °C in the dark, 32 isolates were randomly picked per rhizosphere soil sample. To avoid potential fungal contamination, only highly diluted samples were used for isolation. The isolates were then re-streaked on TSA plates for colony purification. Approximately 5.5% (35 isolates) of the bacterial isolates failed to grow on the TSA plates for unknown reasons when we re-streaked them and were therefore omitted from the dataset. The final collection thus consisted of 605 bacterial isolates derived from 20 rhizosphere soil samples. All purified isolates were cultured in 100 μl tryptone soy broth (TSB, liquid TSA) in 96-well microtiter plates at 30 °C with shaking (rotary shaker at 170 rpm) for 18 h before freezing and storing at −80 °C in 15% glycerol.Strain identificationWe sequenced the full 16 S rRNA gene to taxonomically identify all 605 rhizobacterial isolates. The 16 S rRNA gene was sequenced via Sanger sequencing of PCR products from glycerol stocks by Shaihai Songon Biotechnology Co., Ltd, Shaihai Station. The PCR system (25 µl) was composed of 1 µl of bacterial cells (overnight culture), 12.5 µl mixture, 1 µl of forward (27 F: 5-AGA GTT TGA TCA TGG CTC AG-3) and reverse primer (1492 R: 5-TAC GGT TAC CTT GTT ACG ACT T-3) each [17] and 9.5 µl of sterilized water. PCR was performed by initially denaturizing at 95 °C for 5 min, cycling 30 times with a 30-s denaturizing step at 94 °C, annealing at 58 °C for 30 s, extension at 72 °C for 1 min 30 s, and a final extension at 72 °C for 10 min. The 16 S rRNA gene sequences were identified using NCBI databases and homologous sequence similarity. A total of 90 bacterial isolates that were identified as Ralstonia solanacearum were removed from further analyses, resulting in 515 remaining isolates.Direct effect of rhizobacteria on pathogen growth in vitroWe used R. solanacearum strain QL-Rs1115 tagged with the pYC12-mCherry plasmid as a model bacterial pathogen [8, 18]. We first tested the direct effects of the 515 non-R. solanacearum bacterial strains on the growth of R. solanacearum in vitro by using supernatant assays. Briefly, after 48 h of growth in NB (nutrient broth) medium (glucose 10.0 g l−1, tryptone 5.0 g l−1, yeast extract 0.5 g l−1, beef extract 3.0 g l−1, pH 7.0) on a shaker at 170 rpm, 30 °C, all bacterial cultures were filter sterilized to remove living cells (0.22 µm filter). Subsequently, 20 µl of sterile supernatant from each strain’s culture and 2 µl overnight culture of the pathogen (adjusted to OD600 = 0.5 after 12 h growth at 30 °C with shaking) were added into 180 µl of fresh NB medium (5-times diluted, in order to better reflect the effect of the supernatant). Control treatments were inoculated with 20 µl of 5 X diluted NB media instead of the bacterial supernatant. Each treatment was conducted in triplicate. All bacterial cultures were grown for 48 h at 30 °C with shaking (170 rpm) before measuring pathogen density as red mCherry protein fluorescence intensity (excitation: 587 nm, emission: 610 nm) [9, 11] which was linearly related to the CFU of pathogen R. solanacearum (Fig. S1). To test for significance of growth promotion or inhibition, R. solanacearum densities were log10-transformed prior to analyses of variance (ANOVA) and Bonferroni t test to compare mean differences between each rhizobacterial supernatant treatment and the control treatment, with p values less than 0.05 considered statistically significant. The effect on pathogen growth was defined as the percentage of improvement or reduction in pathogen growth by the supernatant compared to the control treatment. When the effect on pathogen growth was positive, i.e., when the supernatants from strains significantly promoted the growth of the pathogen, they were considered as helpers of the pathogen. If the effect on pathogen growth was negative, i.e., when the supernatants from strains significantly inhibited the growth of the pathogen, they were considered as inhibitors of the pathogen.Assessing strain redundancy among the 515 non-Ralstonia solanacearum bacteriaWe assessed possible redundancy among the 515 strains of the non-Ralstonia solanacearum rhizobacteria. To encompass both taxonomic and functional redundancies, we considered the 16 S rRNA gene sequences as well as the direct effect of their supernatant on Ralstonia solanacearum. Self BLAST searches were performed on the full 515 sequence dataset using the makeblastdb and blastn commands from the BLAST command line tool [19]. Sequences showing >99% identity over >95% of the full length of the 16 S rRNA gene were considered as taxonomically redundant. We then compared the direct effects on pathogen growth of the taxonomically redundant strains, and removed those showing the same patterns of interactions (positive, negative or neutral). Accordingly, (see the dataset “Library of rhizobacterial strains” in the supplementary information), 355 of the 515 strains (68.9%) were removed from the original dataset for further analyses.Phylogenetic tree constructionThe 16 S rRNA gene sequences of the 160 non-redundant bacteria were aligned using MUSCLE [20]. Sequences in the alignment were trimmed at both ends to obtain maximum overlap using the MEGA X software, which was also used to construct taxonomic cladograms [21]. We constructed a maximum-likelihood (ML) tree, using a General Time Reversible (GTR) + G + I model, which yielded the best fit to our data set. Bootstrapping was carried out with 100 replicates retaining gaps. A taxonomic cladogram was created using the EVOLVIEW web tool (https://evolgenius.info//evolview-v2/). To show the relationship between phylogeny and the effects of rhizobacteria on pathogen growth, we added taxonomic status (phylum) of each rhizobacterial strain and its effect on pathogen growth as heatmap rings to the outer circle of the tree separately (Fig. 2B).Fig. 2: Taxonomic characterization of rhizobacterial isolates that inhibited or helped the growth of Ralstonia solanacearum.A Distribution of in vitro effects of 160 rhizobacterial supernatants on R. solanacearum growth. The red vertical line represents no effect on R. solanacearum growth. B Cladogram depicting the phylogenetic relationship among the 160 isolates based on their full-length 16 S rRNA gene sequences. The inner ring depicts the different effect of isolates supernatant on R. solanacearum growth: positive effect (blue), negative effect (red) and no significant effect (gray). The outer ring shows the four phyla to which the isolates belong. C The proportion of rhizobacterial isolates per phylum whose supernatant showed inhibitory, stimulatory or no effect on R. solanacearum growth. The size of the circles represents the number of rhizobacterial isolates in the given phylum. The thickness of lines represents the percentage of rhizobacterial isolates that have the indicated effect on R. solanacearum growth in each phylum.Full size imageEffects of rhizobacteria on pathogen helper strains growth in vitroWe then assessed the potential of different rhizosphere isolates to inhibit helper strains. We first selected two model helper strains (Phyllobacterium ifriqiyense LM1 (Pi) and Microbacterium paraoxydans LM2 (Mp)), which showed strong positive effects on pathogen growth both in co-culture and in supernatant assays (Fig. S2). We defined the effect of rhizobacterial strains on the growth of helpers as the indirect effect on R. solanacearum growth. To study these indirect effects, we first chose a subset of 46 rhizobacterial strains representing a gradient of positive, neutral or negative effect on pathogen growth based on supernatant assays (results in x axis of Figs. 3C and 4A, B, C). We then tested the effects of these 46 rhizobacterial strains on the growth of each of the two helper strains using supernatant assays. Briefly, after 48 h growth in NB media, each of the 46 bacterial monocultures was passed through a 0.22 µm filter to remove living cells. Then 20 µl of sterile supernatant from each strain’s culture and 2 µl overnight culture of Pi or Mp (adjusted to OD600 = 0.5 after 12 h growth at 30 °C with shaking) were added into 180 µl of fresh NB medium (5-times diluted, in order to better reflect the effect of the supernatant). Control treatments were inoculated with 20 µl of 5× diluted NB media instead of a bacterial supernatant. Each treatment was replicated four times. All bacterial cultures were grown for 24 h at 30 °C with shaking (170 rpm) before measuring helper density as optical density (OD600). To test for significance of growth promotion or inhibition, we used analyses of variance (ANOVA) and Bonferroni t test to compare mean differences of helper density between each rhizobacterial supernatant treatment and the control treatment, with p values lower than 0.05 being considered statistically significant. The effect of rhizobacteria on the helpers’ growth (results in y axis of Fig. 3C and x axis of Fig. 4D, E, F) was defined as the percentage of increase or reduction in helper growth by the supernatant compared to the control treatment.Fig. 3: Effect of helper strains on Ralstonia solanacearum growth and plant disease severity.Effects of the two helper strains Phyllobacterium ifriqiyense (Pi) and Microbacterium paraoxydans (Mp) on Ralstonia solanacearum (Rs) growth in vitro (A) and in vivo (B) and on plant disease severity (C). Different letters indicate significant differences based on Tukey post hoc test. Error bars show ±1 SE (n = 3 for in vitro, n = 4 for in vivo). D Effects of 46 rhizobacterial strains on the growth of R. solanacearum and the two model helper strains in vitro. The x-axis shows the direct effect of each rhizobacterial strain on R. solanacearum growth (data from the experiment in which R. solanacearum was grown in the presence of supernatant from each of the 46 rhizobacterial strains—the same data is presented on the x axis of Fig. 4A). The y-axis shows the effect of each rhizobacterial strain on each of the two helper strains (data from the experiment in which each helper was grown in the presence of supernatant from each of the 46 rhizobacterial strains—the same data is presented on the x axis of Fig. 4C). In (C), “−1”, “0” and “1” on the x-axis denote that R. solanacearum growth is completely inhibited, not influenced or increased 2× by supernatant from the rhizobacteria, respectively. Similarly, “−1”, “0” and “1” on the y-axis denote the same growth effects with reference to growth of the helper strains. Black dots indicate results involving interactions with Pi, and red dots indicate results involving interactions with Mp.Full size imageFig. 4: The importance of direct versus indirect effects on Ralstonia solanacearum density and disease severity in the presence of helper strains.In the presence of helper Phyllobacterium ifriqiyense (Pi) or Microbacterium paraoxydans (Mp), respectively, the importance of direct effects on the density of R. solanacearum both (A) in vitro and (B) in vivo, and (C) disease severity (the data on the x axis of (A) are the same data which was presented on the x axis of Fig. 3C, the data on x axis of (B) and (C) are part of the data on x axis of (A)); the importance of indirect effects on the density of R. solanacearum both (D) in vitro and (E) in vivo, and (F) disease severity (the data on the x axis of (D) are the same data which was presented on the y axis of Fig. 3C, the data on x axis of (E) and (F) are part of the data on x axis of (D)). In all panels, “−1”, “0” and “1” on the x-axis denote that R. solanacearum growth (A, B, and C) or helper growth (D, E, and F) is completely inhibited, not influenced or increased 2× by supernatant from the rhizobacteria, respectively.Full size imageIn vitro pathogen growth in the presence of a helper strain and supernatant from rhizobacterial isolatesTo disentangle the direct effects from the indirect effects of rhizobacteria on R. solanacearum growth, we compared their relative effects using in vitro triculture assays comprised of R. solanacearum, one of the two helper strains and supernatant of one of the 46 chosen rhizobacterial strains. Briefly, after 48 h of growth in NB media, each of the 46 bacterial monocultures was passed through a 0.22 µm filter to remove living cells. Then, 20 µl of sterile supernatant from each strain’s culture and 2 µl overnight culture of Pi or Mp (densities were adjusted to ~107 cells per ml) were added to 180 µl of fresh NB medium (5-times diluted). Each treatment was replicated four times. At the same time, 2 µl overnight culture of mCherry-tagged R. solanacearum (density was adjusted to ~106 cells per ml) was added to each treatment in 96-well plates at 30 °C with shaking (170 rpm). After 24-h growth, R. solanacearum density (results in y axis of Fig. 4A, D) was measured as the red mCherry protein fluorescence intensity (excitation: 587 nm, emission: 610 nm) with a SpectraMax M5 plate reader.In vivo pathogen growth and plant disease development in the presence of a helper strain and a rhizobacterial strainTo validate in vitro results, we set up greenhouse experiments where plants were inoculated with a bacterial consortium consisting of R. solanacearum, one of the two helper strains and a test rhizobacterial strain. Tomato seeds (Lycopersicon esculentum, cultivar “Ai hong sheng”) were surface-sterilized by soaking them in 3% NaClO for 5 min and in 70% ethyl alcohol for 1 min before being germinated on water-agar plates for 2 days. Seeds were then sown into seedling trays containing gamma irradiation-sterilized (to avoid potential effects of the resident community) seedling substrate (Huainong, Huaian Soil and Fertilizer Institute). At the three-leaf stage, tomato plants were transplanted to seedling trays containing 200 g of the same seedling substrate as describe above.To relate our results to practical application conditions, we selected a subset of 12 strains that displayed a range of inhibitions effects on pathogen and helpers (Table S1) out of the 46 rhizobacterial isolates used for the in vitro assays. Each rhizobacterial strain was used in combination with each of the two helper strains and R. solanacearum, resulting in a total of 28 treatments (Table S2), including a water control, R. solanacearum alone, and R. solanacearum with just each of the two helper strains (results in Fig. 3B, C). For each treatment, four replicate seedling trays were used, with each replicate seedling tray containing 4 tomato plants. Three days after transplantation, plants of each treatment were inoculated with one of the two helper strains, alone or in combination with one of the rhizobacterial strains, using the root drenching method at a final concentration of 108 CFU g−1 soil for each bacterial strain [22]. Seven days after inoculation of helper alone or together with rhizobacteria, R. solanacearum was introduced to the roots of all plants at a final concentration of 107 CFU g−1 soil. The positive control treatment with R. solanacearum alone was inoculated only with the pathogen, and the negative control treatment was not inoculated with any bacteria. Tomato plants were maintained under standard greenhouse conditions (i.e., at natural temperature variation ranging from 28 °C to 32 °C, 15/9 h day/night conditions) and watered regularly with sterile water. Seedling trays were rearranged randomly every two days. Forty days after transplantation, plants were destructively harvested. The disease index for each plant was recorded based on a scale ranging from 0 to 4 [23]. Disease severity for each replicate seedling plate was calculated as described by: Disease severity = [∑ (The number of diseased plants in the disease index category × disease index category)/ (Total number of plants used in the experiment × highest disease index category)] ×100% [23, 24]. Simultaneously, we collected rhizosphere soil samples following an established protocol [4]. Briefly, two plants were randomly chosen from each replicate seedling tray to collect rhizosphere soils and further combined to yield one sample, resulting in a total of 112 rhizosphere soil samples for which R. solanacearum population densities were determined.Quantification of R. solanacearum at the end of the in vivo experimentWe determined R. solanacearum densities using quantitative PCR (qPCR). DNA was extracted from rhizosphere soils using a Power Soil DNA isolation kit (Mo Bio Laboratories) following the manufacturer’s protocol. DNA concentrations were determined by using a NanoDrop 1000 spectrophotometer (Thermo Scientific) and extracted DNA was used for R. solanacearum density measurements using specific primers (forward, 5ʹ-GAA CGC CAA CGG TGC GAA CT-3ʹ; reverse, 5ʹ-GGC GGC CTT CAG GGA GGT C-3ʹ) targeting the fliC gene, which encodes the R. solanacearum flagellum subunit [25]. The qPCR analyses were carried out with a StepOnePlus Real-Time RCR Instrument using SYBR green fluorescent dye detection and three technical replicates as described previously [4].Statistical analysesTo meet assumptions of normality and homogeneity of variance, R. solanacearum densities measured in vitro and in vivo were log10-transformed. When comparing mean differences between treatments, we used analyses of variance (ANOVA) and the Tukey Test, where p values lower than 0.05 were considered statistically significant. R. solanacearum densities were explained by two quantitative indices, the direct effect of rhizobacteria on R. solanacearum growth (the effect of rhizobacteria on R. solanacearum growth) and the indirect effect of rhizobacteria on R. solanacearum growth (the effect of rhizobacteria on helper strains’ growth). Nonlinear regression analyses (Sigmoidal, Sigmoid, 3 Parameter) were used to analyze the relationship between the direct effect and pathogen density, as well as the relationship between indirect effects and pathogen density in the presence of helper strains in vitro. The relationships between them, and between direct/indirect effects and disease severity in the presence of helper strains in vivo, were analyzed using linear regressions. These analyses were carried out using the R 3.6.3 program (www.r-project.org) and Sigma Plot (V.12.5).To further consider the growth inhibition of R. solanacearum, and disease suppression, we fitted a linear model to estimate the relative importance of direct effects versus indirect effects on the density of R. solanacearum both in vitro and in vivo, and on disease severity. This model considered the interaction scenario where rhizobacterial strains inhibited both the pathogen and its helpers (see the R script “Model” in the supplementary information). These analyses were performed in R version 3.6.3 [26] in conjunction with the package car, readxl and dplyr, and tidyverse 1.2.1 [27]. Briefly, proportional effects were normalized using a folded cube root transformation as suggested in J.W. Tukey [28] and fitted using a linear model with direct effects, indirect effects, and an interaction between helper strains and indirect effects as fixed factors. Normality of residuals was tested using the Shapiro-Wilk normality test and visual inspection of QQ-plots with standardized residuals. Type-II sum of squares were calculated using the ANOVA function from car 3.0-2 [29]. Subsequent visualization of the model outcome (results in Fig. 5) showed the predicted R. solanacearum densities and disease severity for different values of the inhibition via pathogen (Direct) or helper (Indirect) as estimated from the statistical model. For the Direct effect line, the indirect effect is set to be zero, while for the Indirect effect line, the direct effect is set to be zero.Fig. 5: The relative importance of direct versus indirect effects on Ralstonia solanacearum density and disease severity in the presence of helper strains.Relative importance of direct versus indirect effects on Ralstonia solanacearum density both in vitro (A) and in vivo (B), and disease severity (C) in presence of helper strains on the interaction scenario where rhizobacterial strains inhibited both the pathogen and its helpers (quadrant “H−P−” in Fig. 3C). This shows the predicted R. solanacearum densities and disease incidence for different values of the inhibition via pathogen (Direct) or helper (Indirect) as estimated from the statistical model (Table 1) which with direct effects, indirect effects, and an interaction between helper strains and indirect effects as fixed factors. For the Direct line, the indirect effect was set to zero, while for the indirect line, the direct effect was set to zero.Full size image More

NEWS AND VIEWS
19 October 2021

Rhinoceros genomes uncover family secrets

Genomes from living and extinct rhinos reveal that different species evolved as a result of geographic isolation. A comparison of DNA from different species also shows that rhinos have long displayed low genetic variability.

Desire Lee Dalton

0
&

Stefan Prost

1

Desire Lee Dalton

Desire Lee Dalton is at the South African National Biodiversity Institute, Pretoria 0001, South Africa.

View author publications

You can also search for this author in PubMed
 Google Scholar

Stefan Prost

Stefan Prost is at the South African National Biodiversity Institute, Pretoria 0001, South Africa, and in the Department of Behavioural and Cognitive Biology, University of Vienna, the Konrad Lorenz Institute of Ethology at the Vetmeduni Vienna, and the Natural History Museum, Vienna, Austria.

View author publications

You can also search for this author in PubMed
 Google Scholar

Share on Twitter
Share on Twitter

Share on Facebook
Share on Facebook

Share via E-Mail
Share via E-Mail

Historically, rhinos were once abundant throughout Europe, Asia and Africa1. Today, five species of rhinoceros survive as small populations in Asia and Africa, and are all threatened with extinction2. Although well studied, there is debate in the literature about evolutionary relationships between modern and extinct rhinos, with three hypotheses being proposed (Fig. 1a–c). Writing in Cell, Liu et al.3 analyse contemporary and ancient rhinoceros DNA to piece together the puzzle of the rhino’s evolutionary history.

Access options

Access through your institution

Change institution

Buy or subscribe

/* style specs start */
style{display:none!important}.LiveAreaSection-193358632 *{align-content:stretch;align-items:stretch;align-self:auto;animation-delay:0s;animation-direction:normal;animation-duration:0s;animation-fill-mode:none;animation-iteration-count:1;animation-name:none;animation-play-state:running;animation-timing-function:ease;azimuth:center;backface-visibility:visible;background-attachment:scroll;background-blend-mode:normal;background-clip:borderBox;background-color:transparent;background-image:none;background-origin:paddingBox;background-position:0 0;background-repeat:repeat;background-size:auto auto;block-size:auto;border-block-end-color:currentcolor;border-block-end-style:none;border-block-end-width:medium;border-block-start-color:currentcolor;border-block-start-style:none;border-block-start-width:medium;border-bottom-color:currentcolor;border-bottom-left-radius:0;border-bottom-right-radius:0;border-bottom-style:none;border-bottom-width:medium;border-collapse:separate;border-image-outset:0s;border-image-repeat:stretch;border-image-slice:100%;border-image-source:none;border-image-width:1;border-inline-end-color:currentcolor;border-inline-end-style:none;border-inline-end-width:medium;border-inline-start-color:currentcolor;border-inline-start-style:none;border-inline-start-width:medium;border-left-color:currentcolor;border-left-style:none;border-left-width:medium;border-right-color:currentcolor;border-right-style:none;border-right-width:medium;border-spacing:0;border-top-color:currentcolor;border-top-left-radius:0;border-top-right-radius:0;border-top-style:none;border-top-width:medium;bottom:auto;box-decoration-break:slice;box-shadow:none;box-sizing:border-box;break-after:auto;break-before:auto;break-inside:auto;caption-side:top;caret-color:auto;clear:none;clip:auto;clip-path:none;color:initial;column-count:auto;column-fill:balance;column-gap:normal;column-rule-color:currentcolor;column-rule-style:none;column-rule-width:medium;column-span:none;column-width:auto;content:normal;counter-increment:none;counter-reset:none;cursor:auto;display:inline;empty-cells:show;filter:none;flex-basis:auto;flex-direction:row;flex-grow:0;flex-shrink:1;flex-wrap:nowrap;float:none;font-family:initial;font-feature-settings:normal;font-kerning:auto;font-language-override:normal;font-size:medium;font-size-adjust:none;font-stretch:normal;font-style:normal;font-synthesis:weight style;font-variant:normal;font-variant-alternates:normal;font-variant-caps:normal;font-variant-east-asian:normal;font-variant-ligatures:normal;font-variant-numeric:normal;font-variant-position:normal;font-weight:400;grid-auto-columns:auto;grid-auto-flow:row;grid-auto-rows:auto;grid-column-end:auto;grid-column-gap:0;grid-column-start:auto;grid-row-end:auto;grid-row-gap:0;grid-row-start:auto;grid-template-areas:none;grid-template-columns:none;grid-template-rows:none;height:auto;hyphens:manual;image-orientation:0deg;image-rendering:auto;image-resolution:1dppx;ime-mode:auto;inline-size:auto;isolation:auto;justify-content:flexStart;left:auto;letter-spacing:normal;line-break:auto;line-height:normal;list-style-image:none;list-style-position:outside;list-style-type:disc;margin-block-end:0;margin-block-start:0;margin-bottom:0;margin-inline-end:0;margin-inline-start:0;margin-left:0;margin-right:0;margin-top:0;mask-clip:borderBox;mask-composite:add;mask-image:none;mask-mode:matchSource;mask-origin:borderBox;mask-position:0% 0%;mask-repeat:repeat;mask-size:auto;mask-type:luminance;max-height:none;max-width:none;min-block-size:0;min-height:0;min-inline-size:0;min-width:0;mix-blend-mode:normal;object-fit:fill;object-position:50% 50%;offset-block-end:auto;offset-block-start:auto;offset-inline-end:auto;offset-inline-start:auto;opacity:1;order:0;orphans:2;outline-color:initial;outline-offset:0;outline-style:none;outline-width:medium;overflow:visible;overflow-wrap:normal;overflow-x:visible;overflow-y:visible;padding-block-end:0;padding-block-start:0;padding-bottom:0;padding-inline-end:0;padding-inline-start:0;padding-left:0;padding-right:0;padding-top:0;page-break-after:auto;page-break-before:auto;page-break-inside:auto;perspective:none;perspective-origin:50% 50%;pointer-events:auto;position:static;quotes:initial;resize:none;right:auto;ruby-align:spaceAround;ruby-merge:separate;ruby-position:over;scroll-behavior:auto;scroll-snap-coordinate:none;scroll-snap-destination:0 0;scroll-snap-points-x:none;scroll-snap-points-y:none;scroll-snap-type:none;shape-image-threshold:0;shape-margin:0;shape-outside:none;tab-size:8;table-layout:auto;text-align:initial;text-align-last:auto;text-combine-upright:none;text-decoration-color:currentcolor;text-decoration-line:none;text-decoration-style:solid;text-emphasis-color:currentcolor;text-emphasis-position:over right;text-emphasis-style:none;text-indent:0;text-justify:auto;text-orientation:mixed;text-overflow:clip;text-rendering:auto;text-shadow:none;text-transform:none;text-underline-position:auto;top:auto;touch-action:auto;transform:none;transform-box:borderBox;transform-origin:50% 50% 0;transform-style:flat;transition-delay:0s;transition-duration:0s;transition-property:all;transition-timing-function:ease;vertical-align:baseline;visibility:visible;white-space:normal;widows:2;width:auto;will-change:auto;word-break:normal;word-spacing:normal;word-wrap:normal;writing-mode:horizontalTb;z-index:auto;-webkit-appearance:none;-moz-appearance:none;-ms-appearance:none;appearance:none;margin:0}.LiveAreaSection-193358632{width:100%}.LiveAreaSection-193358632 .login-option-buybox{display:block;width:100%;font-size:17px;line-height:30px;color:#222;padding-top:30px;font-family:Harding,Palatino,serif}.LiveAreaSection-193358632 .additional-access-options{display:block;font-weight:700;font-size:17px;line-height:30px;color:#222;font-family:Harding,Palatino,serif}.LiveAreaSection-193358632 .additional-login >li:not(:first-child)::before{transform:translateY(-50%);content:”;height:1rem;position:absolute;top:50%;left:0;border-left:2px solid #999}.LiveAreaSection-193358632 .additional-login >li:not(:first-child){padding-left:10px}.LiveAreaSection-193358632 .additional-login >li{display:inline-block;position:relative;vertical-align:middle;padding-right:10px}.BuyBoxSection-683559780{display:flex;flex-wrap:wrap;flex:1;flex-direction:row-reverse;margin:-30px -15px 0}.BuyBoxSection-683559780 .box-inner{width:100%;height:100%}.BuyBoxSection-683559780 .readcube-buybox{background-color:#f3f3f3;flex-shrink:1;flex-grow:1;flex-basis:255px;background-clip:content-box;padding:0 15px;margin-top:30px}.BuyBoxSection-683559780 .subscribe-buybox{background-color:#f3f3f3;flex-shrink:1;flex-grow:4;flex-basis:300px;background-clip:content-box;padding:0 15px;margin-top:30px}.BuyBoxSection-683559780 .title-readcube{display:block;margin:0;margin-right:20%;margin-left:20%;font-size:24px;line-height:32px;color:#222;padding-top:30px;text-align:center;font-family:Harding,Palatino,serif}.BuyBoxSection-683559780 .title-buybox{display:block;margin:0;margin-right:29%;margin-left:29%;font-size:24px;line-height:32px;color:#222;padding-top:30px;text-align:center;font-family:Harding,Palatino,serif}.BuyBoxSection-683559780 .title-asia-buybox{display:block;margin:0;margin-right:5%;margin-left:5%;font-size:24px;line-height:32px;color:#222;padding-top:30px;text-align:center;font-family:Harding,Palatino,serif}.BuyBoxSection-683559780 .asia-link{color:#069;cursor:pointer;text-decoration:none;font-size:1.05em;font-family:-apple-system,BlinkMacSystemFont,”Segoe UI”,Roboto,Oxygen-Sans,Ubuntu,Cantarell,”Helvetica Neue”,sans-serif;line-height:1.05em6}.BuyBoxSection-683559780 .access-readcube{display:block;margin:0;margin-right:10%;margin-left:10%;font-size:14px;color:#222;padding-top:10px;text-align:center;font-family:-apple-system,BlinkMacSystemFont,”Segoe UI”,Roboto,Oxygen-Sans,Ubuntu,Cantarell,”Helvetica Neue”,sans-serif;line-height:20px}.BuyBoxSection-683559780 .access-asia-buybox{display:block;margin:0;margin-right:5%;margin-left:5%;font-size:14px;color:#222;padding-top:10px;text-align:center;font-family:-apple-system,BlinkMacSystemFont,”Segoe UI”,Roboto,Oxygen-Sans,Ubuntu,Cantarell,”Helvetica Neue”,sans-serif;line-height:20px}.BuyBoxSection-683559780 .access-buybox{display:block;margin:0;margin-right:30%;margin-left:30%;font-size:14px;color:#222;opacity:.8px;padding-top:10px;text-align:center;font-family:-apple-system,BlinkMacSystemFont,”Segoe UI”,Roboto,Oxygen-Sans,Ubuntu,Cantarell,”Helvetica Neue”,sans-serif;line-height:20px}.BuyBoxSection-683559780 .price-buybox{display:block;font-size:30px;color:#222;font-family:-apple-system,BlinkMacSystemFont,”Segoe UI”,Roboto,Oxygen-Sans,Ubuntu,Cantarell,”Helvetica Neue”,sans-serif;padding-top:30px;text-align:center}.BuyBoxSection-683559780 .price-from{font-size:14px;padding-right:10px;color:#222;font-family:-apple-system,BlinkMacSystemFont,”Segoe UI”,Roboto,Oxygen-Sans,Ubuntu,Cantarell,”Helvetica Neue”,sans-serif;line-height:20px}.BuyBoxSection-683559780 .issue-buybox{display:block;font-size:13px;text-align:center;color:#222;font-family:-apple-system,BlinkMacSystemFont,”Segoe UI”,Roboto,Oxygen-Sans,Ubuntu,Cantarell,”Helvetica Neue”,sans-serif;line-height:19px}.BuyBoxSection-683559780 .no-price-buybox{display:block;font-size:13px;line-height:18px;text-align:center;padding-right:10%;padding-left:10%;padding-bottom:20px;padding-top:30px;color:#222;font-family:-apple-system,BlinkMacSystemFont,”Segoe UI”,Roboto,Oxygen-Sans,Ubuntu,Cantarell,”Helvetica Neue”,sans-serif}.BuyBoxSection-683559780 .vat-buybox{display:block;margin-top:5px;margin-right:20%;margin-left:20%;font-size:11px;color:#222;padding-top:10px;padding-bottom:15px;text-align:center;font-family:-apple-system,BlinkMacSystemFont,”Segoe UI”,Roboto,Oxygen-Sans,Ubuntu,Cantarell,”Helvetica Neue”,sans-serif;line-height:17px}.BuyBoxSection-683559780 .button-container{display:block;padding-right:20px;padding-left:20px}.BuyBoxSection-683559780 .button-container >a:hover,.Button-505204839:hover,.Button-1078489254:hover{text-decoration:none}.BuyBoxSection-683559780 .readcube-button{background:#fff;margin-top:30px}.BuyBoxSection-683559780 .button-asia{background:#069;border:1px solid #069;border-radius:0;cursor:pointer;display:block;padding:9px;outline:0;text-align:center;text-decoration:none;min-width:80px;margin-top:75px}.BuyBoxSection-683559780 .button-label-asia,.ButtonLabel-3869432492,.ButtonLabel-3296148077{display:block;color:#fff;font-size:17px;line-height:20px;font-family:-apple-system,BlinkMacSystemFont,”Segoe UI”,Roboto,Oxygen-Sans,Ubuntu,Cantarell,”Helvetica Neue”,sans-serif;text-align:center;text-decoration:none;cursor:pointer}.Button-505204839,.Button-1078489254{background:#069;border:1px solid #069;border-radius:0;cursor:pointer;display:block;padding:9px;outline:0;text-align:center;text-decoration:none;min-width:80px;margin-top:10px}.Button-505204839 .readcube-label,.Button-1078489254 .readcube-label{color:#069}
/* style specs end */Subscribe to JournalGet full journal access for 1 year$199.00only $3.90 per issueSubscribeAll prices are NET prices. VAT will be added later in the checkout.Tax calculation will be finalised during checkout.Rent or Buy articleGet time limited or full article access on ReadCube.from$8.99Rent or BuyAll prices are NET prices.

Additional access options:

Log in

Learn about institutional subscriptions

doi: https://doi.org/10.1038/d41586-021-02777-z

References1.Tissier, J. et al. PLoS One 13, e0193774 (2018).PubMed 
Article 

Google Scholar 
2.Ashley, M. V., Melnick, D. J. & Western, D. Conserv. Biol. 4, 71–77 (1990).Article 

Google Scholar 
3.Liu, S. et al. Cell 184, 4874–4885.e16 (2021).PubMed 
Article 

Google Scholar 
4.Van Couvering, J. A. & Delson, E. J. Vert. Paleontol. 40, e1803340 (2020).Article 

Google Scholar 
5.Tougard, C., Delefosse, T., Hänni, C. & Montgelard, C. Mol. Phylogenet. Evol. 19, 34–44 (2001).PubMed 
Article 

Google Scholar 
6.Von Seth, J. et al. Nature Commun. 12, 2393 (2021).PubMed 
Article 

Google Scholar 
7.van der Valk, T., Díez-Del-Molino, D., Marques-Bonet, T., Guschanski, K. & Dalén, L. Curr. Biol. 29, 165–170.e6 (2019).PubMed 
Article 

Google Scholar 
8.Antoine, P.-O. et al. Zool. J. Linn. Soc. 160, 139–194 (2010).Article 

Google Scholar 
9.Cappellini, E. et al. Nature 574, 103–107 (2019).PubMed 
Article 

Google Scholar 
10.Steiner, C. C. & Ryder, O. A. Zool. J. Linn. Soc. 163, 1289–1303 (2011).Article 

Google Scholar 
11.Antoine, P.-O. et al. Zool. J. Linn. Soc. https://doi.org/10.1093/zoolinnean/zlab009 (2021).Article 

Google Scholar 
12.Welker, F. et al. PeerJ 5, e3033 (2017).PubMed 
Article 

Google Scholar 
13.Margaryan, A. et al. Zool. J. Linn. Soc. 190, 372–383 (2020).Article 

Google Scholar 
Download references

Related Articles

Read the paper: Ancient and modern genomes unravel the evolutionary history of the rhinoceros family

Million-year-old DNA provides a glimpse of mammoth evolution

The changing face of birds from the age of the dinosaurs

See all News & Views

Subjects

Genomics

Conservation biology

Evolution

Palaeontology

Latest on:

Genomics

Convergent somatic mutations in metabolism genes in chronic liver disease
Article 13 OCT 21

A census of cell types in the brain’s motor cortex
News & Views Forum 06 OCT 21

Francis Collins to step down at NIH: scientists assess his legacy
News 06 OCT 21

Evolution

How ancient reptiles were streamlined for flight
Research Highlight 18 OCT 21

Finicky no more: ancient snakes ate their way to success
Research Highlight 14 OCT 21

Genomic reconstruction of the SARS-CoV-2 epidemic in England
Article 14 OCT 21

Jobs

Open Rank, Tenure-Track or Tenured Faculty Position in Bioinformatics

Indiana University – Purdue University Indianapolis (IUPUI)
Indianapolis, IN, United States

Senior Bioinformatician Research Scientist

Brigham and Women’s Hospital (BWH)
Boston, MA, United States

Director of the Pediatric Genomics Program, University of Chicago

The University of Chicago (UChicago)
Chicago, IL, United States

Director of Center for Food Animal HealthDepartment:FAES | Animal Sciences

The Ohio State University (OSU)
Wooster, OH, United States More

1.Cole, M., Lindeque, P., Halsband, C. & Galloway, T. S. Microplastics as contaminants in the marine environment: A review. Mar. Pollut. Bull. 62, 2588–2597. https://doi.org/10.1016/j.marpolbul.2011.09.025,22001295 (2011).CAS 
Article 
PubMed 

Google Scholar 
2.Cauwenberghe, V. L., Vanreusel, A., Mees, J. & Janssen, C. R. Microplastic pollution in deep-sea sediments. Environ. Pollut. 182, 495–499. https://doi.org/10.1016/j.envpol.2013.08.013,24035457 (2013).Article 
PubMed 

Google Scholar 
3.Cole, M. et al. Microplastic ingestion by zooplankton. Environ. Sci. Technol. 47, 6646–6655. https://doi.org/10.1021/es400663f,23692270 (2013).ADS 
CAS 
Article 
PubMed 

Google Scholar 
4.Rodrigues, J. P., Duarte, A. C., Santos-Echeandía, J. & Rocha-Santos, T. Significance of interactions between microplastics and POPs in the marine environment: A critical overview. TrAC Trends Anal. Chem. 111, 252–260. https://doi.org/10.1016/j.trac.2018.11.038 (2019).CAS 
Article 

Google Scholar 
5.Brennecke, D., Duarte, B., Paiva, F., Caçador, I. & Canning-Clode, J. Microplastics as vector for heavy metal contamination from the marine environment. Estuar Coast Shelf Sci. 178, 189–195. https://doi.org/10.1016/j.ecss.2015.12.003 (2016).ADS 
CAS 
Article 

Google Scholar 
6.Martins, I., Rodríguez, Y. & Pham, C. K. Trace elements in microplastics stranded on beaches of remote islands in the NE Atlantic. Mar. Pollut. Bull. 156, 111270. https://doi.org/10.1016/j.marpolbul.2020.111270 (2020).CAS 
Article 
PubMed 

Google Scholar 
7.Woodall, L. C. et al. The deep sea is a major sink for microplastic debris. R. Soc. Open Sci. 1, 140317. https://doi.org/10.1098/rsos.140317 (2014).ADS 
Article 
PubMed 
PubMed Central 

Google Scholar 
8.Alomar, C., Estarellas, F. & Deudero, S. Microplastics in the Mediterranean Sea: Deposition in coastal shallow sediments, spatial variation and preferential grain size. Mar. Environ. Res. 115, 1–10. https://doi.org/10.1016/j.marenvres.2016.01.005,26803229 (2016).CAS 
Article 
PubMed 

Google Scholar 
9.Wright, S. L., Thompson, R. C. & Galloway, T. S. The physical impacts of microplastics on marine organisms: A review. Environ. Pollut. 178, 483–492. https://doi.org/10.1016/j.envpol.2013.02.031,23545014 (2013).CAS 
Article 
PubMed 

Google Scholar 
10.Fendall, L. S. & Sewell, M. A. Contributing to marine pollution by washing your face: Microplastics in facial cleansers. Mar. Pollut. Bull. 58, 1225–1228. https://doi.org/10.1016/j.marpolbul.2009.04.025,19481226 (2009).CAS 
Article 
PubMed 

Google Scholar 
11.Andrady, A. L. Microplastics in the marine environment. Mar. Pollut. Bull. 62, 1596–1605. https://doi.org/10.1016/j.marpolbul.2011.05.030 (2011). CAS 
Article 
PubMed 

Google Scholar 
12.Pakula, C. & Stamminger, R. Electricity and water consumption for laundry washing by washing machine worldwide. Energy Effic. 3, 365–382. https://doi.org/10.1007/s12053-009-9072-8 (2010).Article 

Google Scholar 
13.Belzagui, F., Crespi, M., Álvarez, A., Gutiérrez-Bouzán, C. & Vilaseca, M. Microplastics’ emissions: Microfibers’ detachment from textile garments. Environ. Pollut. 248, 1028–1035. https://doi.org/10.1016/j.envpol.2019.02.059,31091635 (2019).CAS 
Article 
PubMed 

Google Scholar 
14.Ziajahromi, S., Drapper, D., Hornbuckle, A., Rintoul, L. & Leusch, F. D. L. Microplastic pollution in a stormwater floating treatment wetland: Detection of tyre particles in sediment. Sci. Total Environ. 713, 136356. https://doi.org/10.1016/j.scitotenv.2019.136356 (2020).ADS 
CAS 
Article 
PubMed 

Google Scholar 
15.Fadare, O. O. & Okoffo, E. D. Covid-19 face masks: A potential source of microplastic fibers in the environment. Sci. Total Environ. 737, 140279. https://doi.org/10.1016/j.scitotenv.2020.140279 (2020).ADS 
CAS 
Article 
PubMed 
PubMed Central 

Google Scholar 
16.Klein, S., Worch, E. & Knepper, T. P. Occurrence and spatial distribution of microplastics in river shore sediments of the rhine-main area in Germany. Environ. Sci. Technol. 49, 6070–6076. https://doi.org/10.1021/acs.est.5b00492,25901760 (2015).ADS 
CAS 
Article 
PubMed 

Google Scholar 
17.Rummel, C. D., Jahnke, A., Gorokhova, E., Kühnel, D. & Schmitt-Jansen, M. Impacts of biofilm formation on the fate and potential effects of microplastic in the aquatic environment. Environ. Sci. Technol. Lett. 4, 258–267. https://doi.org/10.1021/acs.estlett.7b00164 (2017).CAS 
Article 

Google Scholar 
18.Claessens, M., Meester, S., Landuyt, V. L., Clerck, K. & Janssen, C. R. Occurrence and distribution of microplastics in marine sediments along the Belgian coast. Mar. Pollut. Bull. 62, 2199–2204. https://doi.org/10.1016/j.marpolbul.2011.06.030,21802098 (2011).CAS 
Article 
PubMed 

Google Scholar 
19.Song, Y. K. et al. A comparison of microscopic and spectroscopic identification methods for analysis of microplastics in environmental samples. Mar. Pollut. Bull. 93, 202–209. https://doi.org/10.1016/j.marpolbul.2015.01.015,25682567 (2015).CAS 
Article 
PubMed 

Google Scholar 
20.Lenz, R., Enders, K., Stedmon, C. A., Mackenzie, D. M. A. & Nielsen, T. G. A critical assessment of visual identification of marine microplastic using Raman spectroscopy for analysis improvement. Mar. Pollut. Bull. 100, 82–91. https://doi.org/10.1016/j.marpolbul.2015.09.026,26455785 (2015).CAS 
Article 
PubMed 

Google Scholar 
21.Majewsky, M., Bitter, H., Eiche, E. & Horn, H. Determination of microplastic polyethylene (PE) and polypropylene (PP) in environmental samples using thermal analysis (TGA-DSC). Sci. Total Environ. 568, 507–511. https://doi.org/10.1016/j.scitotenv.2016.06.017,27333470 (2016).ADS 
CAS 
Article 
PubMed 

Google Scholar 
22.Shim, W. J., Song, Y. K., Hong, S. H. & Jang, M. Identification and quantification of microplastics using Nile red staining. Mar. Pollut. Bull. 113, 469–476. https://doi.org/10.1016/j.marpolbul.2016.10.049,28340965 (2016).CAS 
Article 
PubMed 

Google Scholar 
23.Lv, L. et al. A simple method for detecting and quantifying microplastics utilizing fluorescent dyes – safranine T, fluorescein isophosphate, Nile red based on thermal expansion and contraction property. Environ. Pollut. 255, 113283. https://doi.org/10.1016/j.envpol.2019.113283 (2019).CAS 
Article 
PubMed 

Google Scholar 
24.Cauwenberghe, V. L., Devriese, L., Galgani, F., Robbens, J. & Janssen, C. R. Microplastics in sediments: A review of techniques, occurrence and effects. Mar. Environ. Res. 111, 5–17. https://doi.org/10.1016/j.marenvres.2015.06.007 (2015).CAS 
Article 
PubMed 

Google Scholar 
25.Isobe, A. et al. An interlaboratory comparison exercise for the determination of microplastics in standard sample bottles. Mar. Pollut. Bull. 146, 831–837. https://doi.org/10.1016/j.marpolbul.2019.07.033,31426225 (2019).CAS 
Article 
PubMed 

Google Scholar 
26.Uddin, S., Fowler, S. W., Saeed, T., Naji, A. & Al-Jandal, N. Standardized protocols for microplastics determinations in environmental samples from the Gulf and marginal seas. Mar. Pollut. Bull. 158, 111374. https://doi.org/10.1016/j.marpolbul.2020.111374 (2020).CAS 
Article 
PubMed 

Google Scholar 
27.Vermeiren, P., Muñoz, C. & Ikejima, K. Microplastic identification and quantification from organic rich sediments: A validated laboratory protocol. Environ. Pollut. 262, 114298. https://doi.org/10.1016/j.envpol.2020.114298 (2020).CAS 
Article 
PubMed 

Google Scholar 
28.Coppock, R. L., Cole, M., Lindeque, P. K., Queirós, A. M. & Galloway, T. S. A small-scale, portable method for extracting microplastics from marine sediments. Environ. Pollut. 230, 829–837. https://doi.org/10.1016/j.envpol.2017.07.017 (2017).CAS 
Article 
PubMed 

Google Scholar 
29.Ankley, G. T., Di Toro, D. M., Hansen, D. J. & Berry, W. J. Technical basis and proposal for deriving sediment quality criteria for metals. Environ. Toxicol. Chem. 15, 2056–2066. https://doi.org/10.1002/etc.5620151202 (1996).CAS 
Article 

Google Scholar 
30.Japanese Geotechnical Society (JGS) 0131 (JIS A1204). Test method for particle size distribution of soils (2009).31.Japanese Geotechnical Society JGS 0121(JIS A1203). Test method for water content of soils (2009).32.BS ISO 11277: Soil quality. Determination of particle size distribution in mineral soil material. Method by sieving and sedimentation (2009).33.BS ISO 1377:Part 2: Clause 3.2: Determination of the moisture content of soils (1990).34.BS EN 933-2: Tests for geometrical properties of aggregates. Determination of particle size distribution. Test sieves, nominal size of apertures (2020).35.ASTM D6913/D6913M – 17: Standard Test Methods for Particle-Size Distribution (Gradation) of Soils Using Sieve Analysis (2004).36.ASTM D2216 – 19: Standard Test Methods for Laboratory Determination of Water (Moisture) Content of Soil and Rock by Mass (2019).37.Maxwell, S. H., Melinda, K. F. & Matthew, G. Counterstaining to separate Nile red-stained microplastic particles from terrestrial invertebrate biomass. Environ. Sci. Technol. 54, 5580–5588. https://doi.org/10.1021/acs.est.0c00711 (2020).ADS 
CAS 
Article 

Google Scholar 
38.Ehlers, S. M., Maxein, J. & Koop, J. H. E. Low-cost microplastic visualization in feeding experiments using an ultraviolet light-emitting flashlight. Ecol. Res. 35, 265–273. https://doi.org/10.1111/1440-1703.12080 (2020).Article 

Google Scholar 
39.Karakolis, E. G., Nguyen, B., You, J. B., Rochman, C. M. & Sinton, D. Fluorescent dyes for visualizing microplastic particles and fibers in laboratory-based studies. Environ. Sci. Technol. Lett. 6, 334–340. https://doi.org/10.1021/acs.estlett.9b00241 (2019).CAS 
Article 

Google Scholar 
40.Penthala, R. et al. Synthesis of azo and anthraquinone dyes and dyeing of nylon-6,6 in supercritical carbon dioxide. J. CO2 Util. 38, 49–58. https://doi.org/10.1016/j.jcou.2020.01.013 (2020).CAS 
Article 

Google Scholar 
41.Prata, J. C., Costa, J. P., Duarte, A. C. & Rocha-Santos, T. Methods for sampling and detection of microplastics in water and sediment: A critical review. TrAC Trends Anal. Chem. 110, 150–159. https://doi.org/10.1016/j.trac.2018.10.029 (2019).CAS 
Article 

Google Scholar 
42.Hengstmann, E. & Fischer, E. K. Nile red staining in microplastic analysis – proposal for a reliable and fast identification approach for large microplastics. Environ. Monit. Assess. 191, 612. https://doi.org/10.1007/s10661-019-7786-4,31489505 (2019).Article 
PubMed 

Google Scholar 
43.Jung, M. R. et al. Validation of ATR FT-IR to identify polymers of plastic marine debris, including those ingested by marine organisms. Mar. Pollut. Bull. 127, 704–716. https://doi.org/10.1016/j.marpolbul.2017.12.061 (2018).CAS 
Article 
PubMed 

Google Scholar  More

Study areaThe study area (78 000 km2) is located between 68–71°N and 20–26°E, with strong climatic gradients, ranging from wet maritime to relatively dry continental, over tens of kilometers. The landscape of this climatically sensitive high-latitude region has been affected by multiple glaciations in the past. It includes the Scandes Mountains near the Arctic Ocean and low-relief areas to the south and east. The majority of the region (52%) is underlain by sporadic permafrost. Continuous and discontinuous permafrost are limited to the highest mountains of the study area (2% and 7%, respectively)17,26. This large proportion of sporadic, typically warm and shallow permafrost in the study area indicates that ground thermal response to climate warming can be rapid27. Our data do not cover low-relief plateaus of continuous permafrost (similar to northern Siberia and Alaska), where the generally high ice content of soil may lead to different and enhanced LSP responses under climate warming (e.g., ice wedge degradation and surface ponding) with altered AGB feedbacks43,44.LSP observationsThe data consist of 2917 study sites (each 25 m × 25 m) and includes previously combined observations (both in-situ [n = 581] and remote-sensing [n = 2336]) of the active surface features of three cryogenic LSP common in the area: cryoturbation, solifluction, and nivation. These LSP are mainly associated with seasonal freeze–thaw processes. Cryoturbation (i.e., frost churning) is a general term for soil movement caused by differential heave, and it creates typical surface features such as patterned ground, frost boils and hummocky terrain5. Solifluction is the slow mass wasting of surficial deposits through frost creep and permafrost flow, where gravitation causes frost-heaved soil to settle downwards during the summer thaw, creating features of lobes and terraces50. In addition, solifluction also includes gelifluction which is a mass wasting process caused by high porewater pressure in unconsolidated surface debris creating similar lobes and terraces5,50. We use the term nivation to collectively designate various weathering and fluvial processes which are intensified and depicted by the presence of snowbeds (which in general are melting in mid-July – late-August) and nivation hollows28,51. We expect the presence of such a snowbed to be an indication of active nivation processes, since in these environments the year-to-year spatial snow patters are fairly consistent31.The rationale behind LSP sampling is described in previous geomorphic studies which served as a basis for the used protocol52,53. Due to the large study domain, study objectives (focus on distribution of active surface features, not on activity itself) and modeling resolution (50 m × 50 m), we used a visual method to estimate the presence/absence of the mapped LSP. We used high-resolution aerial photography (spatial resolution of 0.25 m−2) and targeted field surveys (GPS accuracy ~5 m; Garmin eTrex personal navigator) to construct the LSP dataset. A binary variable (1 = presence, 0 = absence) was assigned to each LSP to indicate their evident activity (or absence). The activity/absence of the LSP was visually estimated based on the evidence in ground surface, indicated by e.g., frost-heaving, cracking, microtopography (e.g., erosional and depositional forms), soil displacement indicative to a process form (e.g., solifluction lobes, patterned ground), changes in vegetation cover and late-lying snow. Such indicators average the LSP activity over several years. Even small areas with slight indication of activity were considered active processes. However, such a protocol based on a visual assessment is susceptible for incorrect activity classification; solifluction may be active despite having a complete vegetation cover19 and the presence of late-lying snowbed, although being a good indication28, does not necessarily mean that active nivation processes are present.Remotely sensed vegetation indexFor obtaining remotely sensed vegetation index for the study area, we employed a maximum-value compositing approach. We downloaded all available clear sky (less than 80% land cloud cover) Landsat OLI 8 images overlapping the study area from June to September between 2013 and 2017 (total of 1086 scenes) from the United States Geological Survey (USGS) database (http:\earthexplorer.usgs.gov). Images were USGS surface reflectance products, which were preprocessed (georeferencing, projection, and atmospheric corrections) by USGS54. Landsat-8 satellite is the latest addition to the Landsat mission that has provided repeated land surface information globally since the 1970’s and is the most commonly used fine-scale satellite system for vegetation mapping. The native resolution of the Landsat OLI sensor is 30 m for the spectral bands used in the image processing steps of this study.Normalized difference vegetation index (NDVI), a widely used spectral index to estimate the amount of green vegetation, was calculated as55:$$({{{{{rm{rho }}}}}}{{{{{rm{NIR}}}}}}-{{{{{rm{rho }}}}}}{{{{{rm{red}}}}}})/({{{{{rm{rho }}}}}}{{{{{rm{NIR}}}}}}+{{{{{rm{rho }}}}}}{{{{{rm{red}}}}}})$$
(1)
where ρNIR and ρred are the surface reflectance for their respective Landsat bands, 0.851–0.879 (mu)m and 0.636–0.673 (mu)m.USGS provides pixel-based quality assessment bands for all surface reflectance products. These bands were used to mask clouds, snow, water, and other low-quality pixels from the individual NDVI scenes. Additionally, if the NDVI images still had unphysical values over 1 or under -1, these pixels and their surroundings of 100 m radius were excluded. We determined maximum values for each 30 m resolution pixel of the study area individually. After masking cloud, snow, and water from the scenes, obvious scattered erroneous NDVI values remained in some scenes. Therefore, we excluded the values outside the pixel-based 95% percentile prior to maximum composite.The CFmask cloud detection algorithm that is used to generate the quality assessment band has clear difficulties in distinguishing small snow patches from clouds. As such, a large portion of late-lying snow beds were repeatedly and incorrectly classified as clouds. Moreover, the CFmask algorithm creates buffers around the cloud pixels54, hence much information was lost around the snow patches that were incorrectly identified as clouds. After these processing steps, some pixels around the extreme late-lying snow beds had still too low number of NDVI records to provide reliable NDVI values for the maximum composite. To fill these small and scattered gaps in the initial maximum NDVI composite, we selected 74 mostly cloud-free scenes between August and September. For these 74 scenes, we manually digitized cloud masks to exclude cloud-contaminated pixels with high certainty. Moreover, every pixel must have passed the following quality checks to be included in the gap-filling composite: not classified as water in the USGS quality assessment band; normalized difference snow index (NDSI) value less than 0.4, and blue band reflectance less than 0.1 (to exclude snow); reflectance of red band between 0.03 and 0.4 (second check for water and snow, and deepest shadows); NDVI between 0 and 0.4 (lower threshold to exclude snow and water contamination; higher threshold to exclude erroneous values, as very late snowbed habitats always have very limited vegetation cover). Additionally, if the NDVI images had unphysical values over 1 or under -1, these pixels and their surroundings (200 m radius) were excluded. Pixels in the 74 selected images which passed these checks, were then used to create a secondary maximum NDVI composite that was used to fill the gaps in the initial maximum NDVI composite. The secondary composite comprised 0.4% of the pixels in the final composite. Among all 2917 LSP observation sites, 2.9% were located within the gaps in the initial maximum NDVI composite, and thus received their maximum NDVI values from the secondary NDVI composite.In the used Landsat data, the nivation sites were not covered by snow, but instead were associated with generally lower AGB values as nival processes affect the vegetation’s structure and composition (Supplementary Table 1).Above-ground biomass dataAbove-ground biomass (AGB) reference data were collected from two regions, with a total of 433 sites that represent an area of > 4000 km2 (Supplementary Fig. 9). The first dataset (hereafter BM region 1; centering to ca. 69°N, 21°E) was collected between 2008 and 2011, and the second dataset (BM region 2; centering to ca. 70°N, 26.2°E) between 2015 and 2017. Both study regions are representative of an arctic and alpine treeline ecotone and include data from mountain birch forest to barren oroarctic tundra56,57.The BM region 1 dataset consists of 309 field sites (each 10 m × 10 m), which are located around eight different massifs covering a wide range of environmental conditions (Supplementary Figs. 9–10). Sampling was performed in transects to cover various aspects of the slope (i.e., topoclimatic conditions), starting from the foothill of the mountain, and ending at the summit. A plot was systematically established at every 20 m increase in elevation and recorded with a GPS device. Four clip-harvest biomass samples (20 cm × 20 cm) were taken 5 m from the plot center in every cardinal direction. Two samples were used in bare mountaintops (north, south). The clip-harvest samples were dried for 48 h at +65 °C, and dry weight was recorded. The sample biomass values were converted to g m-2 and the average sample value was calculated for each site (Supplementary Fig. 9). The original BM region 1 dataset contains forest and treeline plots, but these were excluded from the final analyses due to an incomparable tree sampling strategy with BM region 2, which could introduce uncertainty into biomass estimates.The BM region 2 data were collected from three different massifs having an elevation range from 120 m to 1064 m (Supplementary Fig. 9). The biomasses were sampled from 102 sites (each 24 m × 24 m in size) that were chosen using a stratified sampling to cover gradients of thermal radiation (potential incoming solar radiation), soil moisture (topographic wetness index, TWI) and vegetation zone (forest, treeline, and alpine zones). Radiation and TWI were calculated from a 10 m digital elevation model (DEM, provided by the National Land Surveys of Finland and Kartverket, the Norwegian mapping authority), and assigned to one of three classes based on observation percentiles (breaks at 20% and 80%) leading to total of 27 strata. Vegetation zones were digitized based on aerial imagery. After the first field survey, 22 sites were added to account for vegetation types that were not sufficiently represented by the GIS-based stratification. Thus, the total sample size of the BM region 2 dataset is 124 AGB sites.The same clip-harvest sample protocol was used as in BM region 1; additional samples were also taken from 12 m in every cardinal direction, thus each site had eight AGB samples (Supplementary Fig. 9). Trees with diameter at breast height (DBH) greater than 20 mm were measured from a 900 m2 circular plot, which corresponds to the size of the NDVI product resolution. Large stems (DBH  > 80 mm in the forest and 40 mm at the treeline) were measured from the whole plot, whereas smaller stems were measured from five subplots. Specifically, the center subplot was 100 m2, and the four subplots located at 8 m to every cardinal direction were each 12.5 m2. For the subplot observations, we used a plot expansion factor (900/150 = 6) to generalize the observations for the whole plot assuming a homogeneous forest structure i.e., each subplot stem represents six trees within the 900 m2 plot. A total of 98% of the measured stems were mountain birch (Betula pubescens ssp. czerepanovii), making it the most abundant species in the area. For predicting stem biomass, we used the average of three allometric equations58,59,60, in order to reduce the uncertainty related to the transferability of an individual allometric model. In addition, Populus tremula (1% of the observations) were found on low-altitude south-facing slopes, and Salix caprea (1%) in moist, nutrient-rich sites. Species-specific models61,62 were used to estimate their respective stem biomasses. Individual pines (Pinus sylvestris) were scattered in the area but were not present in any of the sampled plots.The plots of above-ground tree biomass were converted to g m−2 and added to the mean clip-harvest AGB to obtain the total vascular plant AGB for each site. The BM region 1 and BM region 2 datasets were combined, and the NDVI value was extracted from the site center coordinates.Spatial autocorrelation (SAC) is a common property of any spatial dataset and means that observations are related to one another by the geographical distance63. SAC in the model residuals violates the independence assumption commonly required by statistical models and can lead to inflated hypothesis testing and biased model estimates64. To investigate whether the plot-scale AGB data are spatially autocorrelated, we calculated semivariogram which describes the spatial dependency between the observations as a function of distance between the point pairs65. Semivariogram were calculated as:$${{{{{rm{gamma }}}}}}(h)=frac{1}{2N(h)}mathop{sum }limits_{i=1}^{{N}_{h}}{left(Zleft({s}_{i}right)-Zleft({s}_{i}+hright)right)}^{2}$$
(2)
where N(h) denotes the number of data pairs within distance h, and (Zleft({s}_{i}right)) is an observation (or model residual) in location i. For the calculation, we used R package gstat66 (version 2.0-0). A visual inspection of the semivariogram indicated spatial autocorrelation at short distances (“AGB” in Supplementary Fig. 11). Therefore, for the NDVI-AGB conversion, we used a generalized least squares modeling (GLS, as implemented in R package nlme67 [version 3.1-137]) that can explicitly account for SAC in the data. For the modeling, the AGB values were log(x+0.1) transformed. The GLS, where AGB was modeled as a function of NDVI, were fitted assuming an exponential spatial correlation structure:$$gamma left(hright)={c}_{0}+cleft(1-{e}^{-h/a}right)$$
(3)
where ({c}_{0}) is the difference between the intercept and origin (i.e., the “nugget” parameter in geostatistics), (c) is the amount of variance (i.e., the “sill”) and a represents the distance of spatial dependency (i.e., the “range”). The fitted GLS was as follows:$${{log }}left({AGB}right)=-1.038629+9.725572times {NDVI}$$
(4)
The estimated spatial correlation parameters were c0 = 0.516, c = 0.484 and a = 260.605, indicating that the distance of spatial autocorrelation extends to ca. 261 m. The semivariogram for the model residuals indicated a notable reduction in the amount of spatial autocorrelation compared to the AGB data (Supplementary Fig. 11). The fitted model explained 70.6% of the deviance in the data. When the predicted values were converted back to the response scale, the model explained 60.5% of the deviance. Therefore, for the subsequent analyses we use the above-ground biomass estimated by the model.Environmental predictorsIn addition to LSP, we used climate, topography, and soil predictors to model AGB. Gridded monthly average temperatures and precipitation data (1981–2010; spatial resolution 50 m × 50 m) based on a large collection ( > 950) of Fennoscandian meteorological stations were used in a spatial interpolation scheme17. Three climate predictors—growing degree days (GDD, °C, base temperature 5 °C), mean February air temperature (Tfeb, °C) and water balance (WAB, mm)—were calculated from the gridded climate data. WAB is the difference between total annual precipitation and potential evapotranspiration (PET), which was estimated from the monthly air temperature and precipitation data68:$${PET}=58.93times {T}_{{above}0^circ C}/12$$
(5)
These climatic predictors were selected to represent different aspects of climate that are critical for tundra vegetation: heat requirements, cold tolerance and moisture availability. In addition, two local scale topographic predictors were calculated from a DEM (spatial resolution of 50 m × 50 m, provided by the National Land Survey Institutes of Finland, Norway, and Sweden): topographic wetness index69 (TWI, a proxy for soil moisture) and potential annual direct solar radiation70 (MJ cm-2 a-1). Slope angle was initially considered as a potential predictor for AGB but was later omitted due to the strong correlation with TWI (-0.93, P ≤ 0.001). We also calculated peat cover (%) from a digital land cover classification71. Here, the native resolution of 100 m was resampled at 50 m to match the resolution of the climatic and topographic predictors, using nearest-neighbor interpolation. The binary peat cover variable was transformed to a continuous scale using a spatial mean filter of 3 × 3 pixels52. Finally, the topmost soil layer of a global gridded soil database72 was used to obtain pH data. Again, the original resolution of 250 m was also resampled to 50 m resolution using bilinear interpolation.Our fine-scale data revealed strong environmental gradients over the 78,000 km2 study area (Supplementary Table 1), most of which were only moderately inter-correlated (Spearman’s correlation coefficient  More

Airborne hyperspectral and thermal image acquisitionWe scanned over one million olive and almond trees between 2011 and 2019 with an airborne imaging spectroscopy and thermal imaging facility targeting infected and healthy trees in seven different regions located in Apulia (Italy), Majorca (Balearic Islands, Spain), Alicante, Cordoba and Seville (mainland Spain). In olive groves, over 200,000 and 372,000 trees were imaged from Xf and Vd outbreaks, respectively. In almond groves, we scanned over 132,000 trees from Xf outbreaks in Alicante and Majorca. To evaluate the effects induced by abiotic stress on spectral plant traits, we surveyed over 370,000 healthy trees (outside the outbreak areas) comprising olive and almond species subjected to a wide range of water stress conditions.We surveyed these areas with airborne hyperspectral and thermal cameras on board a manned aircraft flying at 500 m altitude above ground, yielding 40 cm and 60 cm spatial resolution, respectively. We used a hyperspectral camera (VNIR model, Headwall Photonics, Fitchburg, MA, USA) collecting 260 bands in the 400–885 nm region at 1.85 nm/pixel and 12-bit radiometric resolution with a frame rate of 50 Hz. With this spectral configuration, we captured imagery at 6.4 nm full-width at half-maximum (FWHM) bandwidth and obtained an instantaneous field of view (IFOV) of 0.93 mrad and an angular field of view (FOV) of 49.82 deg with an 8 mm focal length lens. The hyperspectral sensor was radiometrically calibrated in the laboratory using an integrating sphere (CSTM-USS-2000C Uniform Source System, LabSphere, North Sutton, NH, USA). At the time of flight, we measured aerosol optical thickness at 550 nm using a Sunphotometer (Microtops II S model 540, Solar LIGHT Co., Philadelphia, PA, USA). We then applied the resulting atmospheric correction of the calibrated radiance imagery with the SMARTS model51 to derive surface reflectance spectra. We carried out ortho-rectification of the hyperspectral imagery (PARGE, ReSe Applications Schläpfer, Wil, Switzerland) with readings acquired by the inertial measuring unit on board the airborne platform (IG500 model, SBG Systems, France). We applied spatial binning through object-based image analysis, thus increasing the signal-to-noise ratio (SNR) of the instrument. Additionally, we conducted spectral binning to reduce the number of spectral bands (260 bands at 1.85 nm resolution). SNR reached >300:1 after binning. We acquired high-resolution tree-crown temperature images with a thermal camera (FLIR SC655, FLIR Systems, USA) at 640 × 480 pixels resolution using a 24.6 mm f/1.0 lens, sensitive to the 7.5–14 μm spectral range and sensitivity below 50 mK.We identified individual trees in the high-resolution hyperspectral and thermal images using object-based crown detection and segmentation methods52. We then calculated the mean hyperspectral radiance, reflectance and temperature for each pure tree crown within every orchard under evaluation. We based our image segmentation methods on Niblack53 and Sauvola and Pietikäinen54, which allowed the isolation of tree crowns from the soil and shadow components. The segmentation of each tree crown was assessed visually to ensure a minimum number of pure vegetation pixels were selected within each tree crown and also spectrally to evaluate the purity of the reflectance extracted from the crown to avoid spectral mixture with soil, shadows and background components24,35.Collection of Xf and Vd biotic stress field dataField assessments of Xf- and Vd-infected trees were carried out from outbreaks affecting olive and almond species in the indicated regions of Italy and Spain between 2011 and 201924,35,52. During these campaigns, we performed quantitative PCR (qPCR)55 for Xf in olive and almond (Alicante), recombinase-polymerase-amplification (RPA) using the AmplifyRP XRT + test (Agdia®, Inc., Elkhart, IN)56 for Xf in almond (Majorca) or conventional PCR57 assays for Vd, as well as visual assessments in individual trees of disease incidence (DI) and disease severity (DS). A sample was considered positive if Ct values were ≤36 and amplification curves were exponential. PCR/qPCR data for model analysis were transformed to 0 and 1, for negative and positive results, respectively, and Ct values were not used in the analysis (see Supplementary Table 2 for the PCR/qPCR primer sequences for Vd and Xf). DS was scored using a 0–4 rating scale according to the percentage of the tree crown showing disease symptoms.In Apulia, the Xf-olive database comprised a total of 15 olive groves surveyed during the June 2016 and July 2017 campaigns. Visual assessments for infection were conducted on 7296 trees (3324 in 2016 and 3972 in 2017). In 2016, 1886 symptomatic (and 1438 asymptomatic) trees were surveyed (762 trees labelled as DS = 1; 802 DS = 2; 250 DS = 3 and 72 DS = 4). In 2017, 1365 were reported as symptomatic (and 2607 asymptomatic) (686 DS = 1; 542 DS = 2; 122 DS = 3 and 15 DS = 4). qPCR assays were carried out to diagnose Xf infection in 77 olive trees, whereby 39 trees tested negative (qPCR = 0) and 38 tested positive (qPCR = 1).On the island of Majorca and at the Alicante province, the field-based Xf-almond database comprised a total of 19 almond groves surveyed in 2018 and 2019, respectively. In Alicante, the field surveys covered 83 ha with 9 almond groves consisting of 943 almond trees. During the field campaigns, almond trees were visually assessed to evaluate Xf-induced DI and DS indices. From this analysis, we identified 593 symptomatic trees and 350 asymptomatic trees. Out of all symptomatic trees, 163 were rated as DS = 1, 214 DS = 2, 157 DS = 3, and 59 DS = 4. Furthermore, qPCR analysis was carried out on 226 almond trees to diagnose Xf infection, resulting in 48 non-infected (qPCR = 0) almond trees and 178 infected trees (qPCR = 1). In Majorca, field surveys in July 2019 covered a total of 2803 ha and comprised 10 almond groves. During the field campaigns, visual observations were carried out on over 4048 almond trees to assess DI and DS, yielding 1387 symptomatic and 2661 asymptomatic trees. From symptomatic trees, 537 were rated as DS = 1449 DS = 2, 359 DS = 3 and 42 DS = 4. We conducted AmplifyRP XRT + assays on 265 almond trees for diagnosing Xf infection the same day they were sampled and identified 141 negative trees (qPCR = 0) and 124 positive trees (qPCR = 1).We carried out physiological measurements of leaf chlorophyll, anthocyanins, flavonoids and nitrogen contents with a Dualex Scientific + (Force-A, Orsay, France) instrument as well as leaf reflectance (400–1000 nm spectral range) and steady-state chlorophyll fluorescence (Ft) using the PolyPen RP400 and FluorPen FP100 instruments (Photon Systems Instruments, Drasov, Czech Republic) during the field evaluations of almond and olive groves conducted in Majorca, Alicante and Apulia regions. In the Xf-olive study site in Apulia, we generated 1023 leaf measurements with Dualex, 1543 single leaf reflectance spectra, as well as 1402 Ft readings over 67 olive trees. In the Xf-almond study sites in Majorca, we measured 1242 leaves with Dualex, 1094 leaves with the PolyPen and 1218 with the Fluorpen instruments from 87 almond trees across a wide range of disease severity levels. For the Xf-almond study sites located at Alicante, we conducted 1649 leaf measurements with Dualex, 632 leaf measurements with PolyPen and 563 leaf measurements with FluorPen FP100 over 43 almond trees.We assessed Vd-infected olive trees from 11 olive groves by surveying an area of over 3000 ha in Castro del Rio and Ecija, southern Spain, in 2011 and 2013, respectively. In Castro del Rio, we conducted visual assessments in an infected area of 96 ha comprising 1878 olive trees, thus identifying 1569 asymptomatic and 283 symptomatic olive trees. Out of the 283 symptomatic trees, 218 were rated as DS = 1; 45 DS = 2; 12 DS = 3 and 8 DS = 4. We measured leaf Fs and Fm’ fluorescence parameters from 25 leaves per tree using a PAM-2100 Pulse-Amplitude Modulated Fluorometer (Heinz Walz GMBH, Effeltrich, Germany). In addition, leaf PRI570 was measured from 25 leaves per tree using a custom-made PlantPen device (Photon System Instrument, Drasov, Czech Republic). Finally, we measured leaf conductance (Gs) on five leaves per tree using a leaf porometer (model SC-1, Decagon Devices, Washington, DC, USA). In the Écija region, the surveyed area covered 3424 ha, and 5223 olive trees were evaluated. We performed visual assessment to determine DI and DS indices of Vd-infected trees, identifying 5040 asymptomatic olive trees. Of the remaining 183 olive trees that were symptomatic, 112 were trees rated as DS = 1; 41 DS = 2; 22 DS = 3 and 8 DS = 4.Trees were evaluated for disease severity and incidence by visual assessment in each outbreak region. PCR assays were carried out on a subset of these trees within each orchard to (i) validate that the pathogen (Xf or Vd) was actually present and the biotic source of symptoms; and (ii) validate that asymptomatic (DS = 0) but infected (PCR = 1) trees were detected using the hyperspectral plant traits estimated through the methodology described in this paper. In general, PCR assays are (i) time consuming and costly, and (ii) difficult to make in large infected trees due to the non-uniform distribution of the infection within each tree crown. These PCR data for each tree along with the field evaluations of DS, DI and non-destructive physiological measurements derived for each tree within every orchard were matched with the high-resolution hyperspectral images to build the biotic databases used in this study. We carried out the field work at each orchard guiding the evaluations and measurements using a high-resolution image to map the location of each tree within the orchard. Due to the planting grids typical of almond and olive species, which were not contiguous or in row-structured patterns, the identification of each individual tree in the images was straightforward.Collection of abiotic stress field dataWe monitored over 3600 ha of olive and almond groves located outside any infected area in Cordoba and Seville, Southern Spain. We performed a multitemporal analysis to study the spectral plant-trait alterations induced by abiotic stress relative to healthy olive and almond trees with data we collected over a 468 ha area comprising two olive and one almond groves throughout July 2016 and August 2017 growing seasons. We analysed 2975 olive and 1964 almond trees in 2016, and 2865 olive and 2063 and almond trees in 2017. At both study sites, we monitored the midday stem water potential (SWP) using a pressure chamber (Soil Moisture Equipment Corp. model 3000, Santa Barbara, CA, USA) on 16 trees per grove. SWP values showed differences between two existing irrigation levels (well-watered and mild water stress), averaging –1.7 and –1.9 MPa across the season in the case of almonds. In olive, SWP in one of the groves reached –3.8 and –3.5 MPa. In 2017, water potential levels averaged –2.9 and –2.0 MPa. In the second grove, irrigation levels were higher, reaching an average SWP of –1.5 MPa. We used an additional study site located in Casariche (Seville province), southern Spain, to validate the results obtained from the multitemporal analysis. This study site covered 3371 ha containing 55 olive groves grown under irrigated and rainfed conditions, with 21,071 olive trees used for statistical analysis.The multitemporal dataset was used to evaluate the water-induced abiotic stress by quantifying the evolution of the importance of the most sensitive spectral traits by clustering non-stressed trees (C0) against groups of trees exposed to increasing levels of water stress (C1 to C4). The multitemporal component of this assessment enabled the evaluation of every single tree across time, therefore selecting the trees for each cluster based on a sustained water stress level, avoiding the selection of trees under short-term stress dynamics. Thus, the clusters were determined based on their CWSI levels, and only the trees with stable water stress levels across two consecutive years (2016 and 2017) were selected for the analysis. For this purpose, we did not include trees that deviated beyond 95% of the CWSI differences calculated between the first and second year in the analysis. After this trimming step, we retained 5484 olive trees (from 5566 trees) and 3652 almond trees (from 3882 almond trees). Trees were then grouped through CWSI clustering analysis using a modified three-sigma rule58. This rule describes the density of a distribution within standard deviation bands on both sides of the mean point into the 68th, 95th and 99.7th percentiles58, representing µ ± σ, µ ± 2σ and µ ± 3σ, respectively. The first interval defined by the classic three-sigma rule (µ ± σ) represented most trees, while the third interval (µ ± 3σ) consisted of very few trees, raising issues for the determination of statistical significance analysis. Based on this observation, we adjusted the breakpoints between groups as follows: we classified those trees that were in the lowest 10th percentile as C0. Trees between the 10th and 68th percentiles (µ + σ) were classified as C1, trees between the 68th and 85th percentile were classified as C2, trees between the 85th and 95th percentile were classified as C3 and finally the trees over the 95th (µ + 2σ) percentile were classified as C4. We thus selected 488 C0, 3066 C1, 1090 C2, 618 C3 and 222 C4 olives trees. Likewise, we grouped almond trees into 390 C0, 1776 C1, 1248 C2, 214 C3 and 24 C4 clusters. Moreover, the analysis of the contribution of a given trait was performed using ML modelling strategies to classify unstressed trees against the clusters defined above that were exposed to increasing levels of water stress. Furthermore, we assessed the consistency of the obtained indicators by performing the classification between stressed and non-stressed trees at an independent olive study site. For this purpose, we evaluated our predictors and compared their contribution over an additional site (Casariche).Model inversion methods for plant-trait estimationWe quantified chlorophyll content (Ca+b), carotenoid content (Cx+c), anthocyanin content (Anth.), mesophyll structure (N), leaf area index (LAI) and average leaf angle (leaf inclination distribution function or LIDF) by radiative transfer model inversion of PROSPECT-D59 and 4SAIL60, as in Zarco-Tejada et al.24. We inverted PROSPECT-D + 4SAIL using a look-up-table (LUT) generated with randomised input parameters. The LUT was generated with 100,000 simulations within fixed ranges (Supplementary Table 3). We implemented a wavelet analysis61 into six wavelets by a Gaussian kernel, estimating the parameters in the top 1% entries ranking the lowest root mean square error (RMSE) values. We then retrieved each plant trait independently by training supported vector machine (SVM) algorithms using the simulated reflectance data as input. We built SVMs in Matlab (MATLAB; Statistics and Machine Learning toolbox and Deep Learning toolbox; Mathworks Inc., Matick, MA, USA) using a Gaussian kernel (radial basis function) with hyperparameters optimised for each model. The training processes were carried out in parallel using the Matlab parallel computing toolbox. With these trained models, we then used the spectral reflectance extracted from the delineated crowns (as show in Fig. 1) to predict plant traits for each individual tree at each study site. The model inversions were carried out for each tree using the crown reflectance. The latter was calculated as an average across all the pixels belonging to the tree crown, delineated using segmentation. This method52 avoids the problem of pixels from within-crown shadows, from tree edges or from sunlit or shaded soil background affecting the spectra, as it retrieves the plant traits from pure sunlit vegetation components of the trees. We also calculated narrow-band spectral indices from reflectance spectra (Supplementary Table 1), which are sensitive to leaf traits and potentially related to disease-induced symptoms. Tree-crown radiance and temperature were used to calculate sun-induced chlorophyll fluorescence at 760 nm (SIF@760) and the crop water stress index (CWSI)37. SIF@760 was quantified using the O2-A in-filling Fraunhofer Line Depth (FLD) method63 and CWSI was calculated by incorporating the tree temperature and the weather data obtained at each study site37.Statistical analysisWe implemented random forest (RF)64 algorithms to classify healthy vs. infected (biotically stressed) trees, and non-stressed vs. water (i.e. abiotically) stressed trees for both tree species. RF algorithms have been widely used in remote sensing studies since they have shown excellent classification accuracies and high processing speeds with high-dimensional data62 and have shown to be accurate in detection of several diseases29,65,66,67. The spectral plant traits estimated by radiative transfer model inversion (Ca+b, Cx+c, Anth., LAI and LIDF), CWSI and SIF@760 were used as inputs for the models. In addition, using a recursive feature elimination approach68 the narrow-band indices that improved the classification in terms of overall accuracy (OA) and kappa coefficient (κ) were added to the models. The pool of narrow-band indices was reduced based on a variance inflation factor (VIF) analysis69 to avoid collinearity among the input features.The RF algorithms were built in Matlab and the hyperparameters were optimised using Bayesian optimisation. The importance of a feature using the RF algorithm was assessed based on the permutation of out-of-bag (OOB) predictor methodology70. To compare the relative differences of the spectral traits in classification of the biotic and abiotic stress, the importance was normalised by dividing the importance of each trait by the highest contribution obtained for each pathogen/species. For the RF models, 500 iterations were run by randomly partitioning each dataset into training (80% of samples) and testing sets (20% of samples). For the training subset, a balanced number of samples from each class was randomly selected at each iteration. The importance obtained by the OOB permutation algorithms was used to build a feature-weighted random forest algorithm (based on Liu and Zhao45), accounting for the importance of each variable on the classification process, evaluating the model against PCR data and visual observations for each biotic stress dataset in terms of OA and κ levels.Probabilities of the predictions were obtained for each sample71 and the uncertain trees were assessed. To extract the uncertainty for each individual tree on the classification, we evaluated the probability distribution for each class from each dataset independently. Then, those trees with a classification probability below the 68th percentile (µ [mean] + σ [standard deviation]) were considered as uncertain and incorporated into a second-stage classification process. The second stage consisted of an unsupervised graph theory–based spectral clustering algorithm72 and included traits selected by focusing on the divergent biotic–abiotic stress obtained from the biotic and the abiotic stress databases. Spectral clustering was performed in R using the kernlab package73.To determine the spectral traits that differed between Xf- and Vd-infected plants and those from the abiotic pathway, we first normalised the importance of the specific traits independently. Then, we compared the common traits between abiotic and biotic stress sets, selecting only biotic stress-related traits that differed in ratio by >0.5 over their homologous abiotic stress trait values. Traits that were only expressed under biotic stress conditions and that showed a normalised importance over 0.5 were included for the second-stage classification process only including those divergent-specific biotic and abiotic stress-related spectral traits as inputs. Specifically, NPQI, Anth. and SIF@760 were considered for the classification of Xf-infected olive trees. Ca+b, SIF@760 and PRIn were used for classifying Xf-infected almond trees. Furthermore, NPQI, Anth. and B spectral traits were selected for classifying uncertain Vd-infected olive trees. Finally, we validated our feature-weighted methodology coupled with the second-stage spectral clustering method against qPCR assays and visual assessment of symptom severity.Reporting summaryFurther information on research design is available in the Nature Research Reporting Summary linked to this article. More

Trends in global biomass C stocksThe CRAFT model reliably reproduces the observed trends in primary and managed forest biomass C stocks (including both above-ground and belowground biomass) in 1990–2020 with a relative root mean square error (RMSE) of 0.57% between simulated and observed biomass C stocks by the FRA2 at the global level. These low divergences between stock estimates result, however, in global C emissions c. 2 times lower according to the CRAFT simulations than the estimates derived from the FRA (Supplementary Table 2). Still, the CRAFT simulations corroborated the FRA observations while adding information on annual estimations of forest C stocks, rather than the 5-years interval data provided by the FRA (Fig. 1a), and dynamic annual net C emissions (Fig. 1c, d) from managed and primary forests. The five sensitivity analyses carried out on the most uncertain model inputs and assumptions (see ‘Methods’ descriptions and Supplementary Figs. 5–10) confirmed the results presented in Fig.1: The largest deviation derived from the sensitivity analysis considering forest gross instead of net area changes results in a relative RMSE of 1.94% with global C emissions c.2.5 times higher than the FRA estimates (Supplementary Table 2). The simulations from the reference model assumptions yield the best RMSE and closest agreement with the C budgets derived from the FRA, indicating that they are the most optimal.Fig. 1: Global trends in total, primary, and managed forests.a Forest biomass C stocks (GtC); b cumulated change in forest area (Mha; negative values indicate area loss); c cumulated C net emissions (GtC; positive values indicate a C source while negative values indicate a C sink); and d cumulated net change in C-stock densities (tC/ha). See Supplementary Fig. 1 for annual fluxes.Full size imageIn line with the FRA data, we find here that the main trend is a loss of total biomass C stocks following three phases: increase in annual emissions, stagnation and slight recovery of C stocks, resulting in net C emissions from forest biomass (Fig. 1c) by 0.74 GtC or 0.03 GtC/yr between 1990 and 2020, contrasted by an opposite trend of increasing biomass density from 70 to 73 tC/ha in total forest (Fig. 1b, d). These figures are within the range of the estimated sink in forest soil and biomass of 0.1 ± 7.3 GtC/yr in 2001–2019 found by Harris et al.17. Our estimation is also consistent with that of Tubiello et al.3 of 0.11 GtC/yr net C emissions from forest ecosystems. A comparison3 of FRA-derived global forest C emissions with other independent estimates reported in 1990–2015 by National Greenhouse Gas Inventories (NGHGIs)—including the Russian Federation, the USA, China, Indonesia, and India—and by the United Nations Framework Convention on Climate Change for other countries (UNFCCC, 202018) yields a slight difference of c. 18%, although the UNFCCC and NGHGI’s account, by definition, only for emissions from managed land3. Further independent comparisons at the national and macro-regional levels are compiled in Supplementary Table 1 and reveal that C emissions estimated in the present study are in good agreement with other research.Here we find that the net C emissions mostly arise from primary forests, which undergo area loss, but also biomass thickening (Fig. 1b, d). By contrast, in spite of area loss, managed forests act as C-sinks following biomass thickening (Fig. 1b, d). Increasing biomass density is therefore key to counteract net C emissions from forest biomass in 1990–2020. While both harvest rate and burnt area increase globally over the period of observation, the increased forest growth rate that we calculate with CRAFT for both primary and managed forests over 1990–2020 emerges here as the only factor explaining increased biomass density at the global level. This is in line with other research pointing to the relevance of biomass thickening for forest C sequestration19. In addition, our finding that the forest growth rate increased annually by 0.19%, 0.21%, and 0.21% from 1990 to 2020, respectively, for primary, managed and total forests of the world is consistent with Kolby Smith et al.20 who find that also net primary production (NPP) increased annually between 0.10 and 0.25% in the period 1982–2011, as well as with other modeling and remote-sensing studies documenting a global greening trend, i.e., vegetation thickening following increased vegetation growth rate21,22. Note that estimates of annual growth rate increase in 1990–2020 by the sensitivity analyses provide narrow ranges of 0.17–0.19, 0.21–0.23, and 0.20–0.22%, respectively, for primary, managed, and total forests of the world (Supplementary Table 2).Proximate drivers of net C emissionsWe develop six counterfactual scenarios23,24,25 in order to investigate how forest biomass density and forest biomass C stocks would evolve in the hypothetical absence of (i) changes in harvest (CF1); (ii) changes in forest growth rates (CF2); (iii) change in burnt area (CF3); (iv) change in forest area (CF4); (v) harvest (CF5); (vi) burnt area (CF6) (see “Methods” section). The comparison of observed and simulated counterfactual trends allows us to isolate and quantify the influence of these four main drivers on global forest C-stock changes at national resolution (CF1 to 4) as well as to quantify the overall effects of total wood extraction and burnt area (CF5 and 6).At the global level, we find that loss of forest area (CF4) is the main driver of the net C emissions from forest biomass (Fig. 2a). In the absence of changes in area, global forest biomass would act as a cumulative net C sink of c. 26.9 GtC in the study period, creating a difference of 27.6 GtC between the actual and the CF4 C budget. This effect in the absence of area change, however, is a composite of an additional C sink of 30.7 in deforesting countries and an additional C source of 3.8 GtC in reforesting countries. Changes in harvest and burnt area from 1990 to 2020 also drove net C emissions from global forest biomass as emissions drop by c. 5.7 and 1.4 GtC in the respective counterfactual scenarios, thus generating net C-sinks of c. 4.9 and 0.63 GtC (Fig. 2a). These figures are in stark contrast with the estimated total sink of c. 49.1 and 5.4 GtC that would emerge in the hypothetical absence of harvest (CF5) and burnt area (CF6; Fig. 2a), respectively. Only changes in forest growth rates counteract the net C emissions from global forest biomass (CF2; Fig. 2a). In the absence of changes in forest growth rates, global forests would act as net C source of c. 7.4 GtC in 1990–2020, i.e., c. 10 times the actually observed source. This net effect in the absence of growth rate change results from an additional C source of 30.4 in countries experiencing growth rate increase and an additional C sink of 23.0 GtC in countries experiencing growth rate decline.Fig. 2: Counterfactual scenarios (1990–2020) assessing the cumulative impact of: changes in harvest (CF1); changes in forest growth rate (CF2); changes in burnt area (CF3); changes in forest area (CF4); total harvest (CF5); and total fire (CF6) on C-dynamics.Panels (a) and (b) show the global country-level gross and net CF C budgets (GtC) and changes in biomass density (tC/ha), respectively, with negative (red) and positive values (blue) indicating net emissions and sinks, respectively, error bars indicate the range of C budgets estimated across the five sensitivity analyses performed to test the model robustness (see Supplementary Fig. 5 for additional figures showing the net difference between CF and actual C budgets and changes in biomass density, Supplementary Table 3 and Supplementary Fig. 5 for results from sensitivity analyses). Maps show the effects of c CF1; d CF2; e CF3; f CF4; g CF5; h CF6, and are represented as the % of actual biomass C stocks that would be reached in each CF in 2020. Values above 100% (red) indicate that actual change result in net C emissions while values below 100% (blue) indicate that actual change result in a net C sink.Full size imageA sensitivity analysis on the potential underestimation of C-dynamics resulting from the use of net area change data at country level (see “Methods” section and Supplementary Fig. 5) reveals that accounting for gross area changes26 instead of net area change would result in higher global C emissions estimates (4.19 GtC in the sensitivity test versus 0.74 GtC in the reference simulation) but would reveal the same patterns of forest C-dynamic drivers (Supplementary Fig. 5). However, the magnitude of the main drivers would be slightly changed with a lower effect of changes in area (C sink in the hypothetical absence of area changes reaching 20.8 GtC in the sensitivity tests versus 26.9 GtC in the reference assessment) and a higher effect of growth rate changes (C source in the hypothetical absence of growth rate changes reaching 13.1 GtC in the sensitivity tests versus 7.4 GtC in the reference assessment). Generally, the range of results derived from the five sensitivity analyses does not change the relative importance of the individual drivers in any of the scenarios (Fig. 2a, Supplementary Table 3, and Supplementary Fig. 5). However, the sensitivity analyses highlight that the uncertainty is large enough to reverse the cumulated C signal in the absence of changes in harvest (CF1), changes in burnt area (CF3), and the complete absence of burnt areas (CF6). By contrast, the signals of CF2 (no growth rate change), CF4 (no area change), and CF5 (no harvest) are larger than the uncertainty across sensitivity analyses, signaling that our findings on these drivers are most robust.The global trends displayed in Fig. 2a, b are the combined results of diverging national forest dynamics (Fig. 2c, h). In particular, shifts in forest area (CF4) contribute to global net C emissions only in the Global South, excluding Vietnam, India, and Chile (Fig. 2f). The impacts of changes in burnt area and harvest are similarly heterogenous, with considerable effects only in some regions (e.g., Vietnam, Mozambique, Fig. 2c, e). In contrast, changes in forest growth rates are more ubiquitous, mainly positive (leading to C-sinks) for most countries, with a few notable exceptions, mainly in arid or boreal regions (e.g., India, Spain, Argentina, Canada; Fig. 2d). Possible reasons explaining the negative effect of change in forest growth rate are forest degradation, increasing drought, cloudiness, or insect outbreaks15,16,17,18,19. Over the period 1990–2020, the strongest harvest impacts are observed in countries with large area of managed forest and high harvest pressure, mostly located in temperate and subtropical areas (CF5; Fig. 2g), while fire impacts are strong in only a few countries (CF6; Fig. 2h).The fact that we use here country-level data comes both with limitations and advantages. The main limitation associated with national data is that it conceals gross C fluxes in forest biomass dynamics and blurs heterogeneity in growth conditions and anthropogenic management within countries. The country-level resolution aggregates the effects of manifold, partly counteracting processes at the local level—including photosynthesis, maintenance respiration, growth respiration, as well as forest area loss and expansion—on the annual dynamic of primary and managed forest biomass. As a consequence, our optimization of the growth function actually reflects apparent national growth rates resulting from the aggregate of these processes. However, this simplification of forest ecosystem functioning is also an advantage. Our approach reproduces forest biomass dynamics very accurately, which is complementary to most process-based models aimed at depicting biological processes and their abiotic controls27 but providing a wide range of C flux estimations1 and hardly reproducing observation from inventory data1,28,29. By contrast, the strength of the modeling approach implemented here is that it can be run with parsimonious data availability and allows to disentangle the major drivers behind forest C-stock and flux trajectories.Typology of forest biomass changeIn order to identify spatial and temporal patterns of drivers in forest biomass trends, we establish a typology of the main drivers over the period 1990–2020 (Fig. 3b). The typology we established is based on the positive versus negative shift in biomass C stocks, and highlights the most important driver of this shift as assessed through the counterfactual assessment, irrespective of the relative importance of the other drivers shown in Fig. 2. However, as the early separation between increasing and decreasing biomass C stocks in the decision tree (Fig. 3b) may conceal the effect of a major driver counteracting the observed C dynamic, the typology also accounts for possible antagonistic effects by identifying cases in which the main driver of observed C-dynamics is not, in absolute terms, the most important driver (e.g., C stocks increase but the driver with the strongest absolute effect counteracts this positive budget, see also Supplementary Fig. 3). By pinpointing the major drivers of forest change at national levels, such an approach enables to identify major levers for forest conservation.Fig. 3: Main drivers of the net C emissions from forest biomass.a Applied at the national level to the 1990–2020 period; b established according to a Boolean typology using the results from the counterfactual scenario assessment as criteria; c enabling to calculate the sum of net C-sinks and net C sources in each type of forest C-dynamics trajectory identified through the typology, error bars indicating the range of C-sinks and sources by main driver estimated across the five sensitivity analyses, with black and gray bars standing, respectively, for solid and hatched countries (see Supplementary Figs. 6–7 for results from sensitivity analyses). The hatches on the countries (a), typology (b), and bar chart (c) stand for cases in which the driver with the strongest effect actually counteracts the observed carbon budget. The color of the hatches corresponds to the main factor identified by the decision tree algorithm. Abbreviation on the typology: E: C sink driven by forest area Expansion; LH: C sink driven by Lower Harvest; FR: C sink driven by Fire Reduction; EG: C sink driven by Enhanced Growth rate; DG: C source driven by Declining Growth rate; FI: C source driven by Fire Increase; HH: C source driven by Higher Harvest; D: C source driven by Deforestation; NS: non-significant change.Full size imageDeforestation was the dominant driver of net C emissions from forest biomass in most countries of South America and Sub-Saharan Africa, corroborating findings from the literature11,30,31 (Fig. 3a, c). The net C emissions by countries where deforestation is the most significant driver reach c. 21.3 GtC, with only 0.3 GtC of these emissions being counteracted by another major driver (either increased growth rate or lower harvest pressure). These emissions represent c. 92.7% of the 21.9 GtC net emissions arising from all countries acting as net C sources (Fig. 3c). Changes in forest growth rates act as the primary drivers in most countries experiencing a net C sink over the period (Fig. 3a, c). The net C-sinks by countries where changes in forest growth rates are the main driver reach c. 16.4 GtC, with 0.9 GtC of these sinks being counteracted by another major driver (increased harvest pressure in all cases except for Sudan where area loss was the major driver counteracting the C sink). These C-sinks mainly driven by increased growth rate represented c. 77.5% of the 21.1 GtC net sink created by all countries acting as net C-sinks (Fig. 3c).Forest area expansion from 1990 to 2020 is the main driver of forest biomass net C sink in only a few Northern countries but also some Southern countries, namely Vietnam, India, and Chile, in line with findings reported for these countries32,33,34, all together accounting for a net C sink of 3.9 GtC. However, more than half of the C-sinks mainly driven by reforestation are counteracted by another major driver (either declining forest biomass growth rate or increased harvest pressure). Similarly, changes in harvest as well as changes in burnt areas are the main drivers of net C sink or source for a handful of countries in 1990–2020 (Fig. 3a). Finally, declining forest biomass growth rate is the primary driver of net C emissions only in Mongolia and Canada, which is consistent with other studies highlighting slower growth, higher mortality, and insect outbreak events in Canadian forests35,36,37.These highlights derived from the typology remained the same in all sensitivity analyses (Supplementary Figs. 6–7), despite some possible changes in country type identification (Fig. 3a and Supplementary Fig. S6) and amplitude shifts in the attribution of main drivers globally (Fig. 3c and Supplementary Fig. 7). The ranges of values in the attribution of main drivers result from the previously reported differences between the counterfactual and actual C budget estimates across sensitivity analyses (see also Supplementary Tables 2–3) combined with some changes in the type of forest C-dynamics trajectory identified through the typology in countries with large forest biomass stocks: China, India, and Australia (Supplementary Note 1 and Supplementary Fig. 6). However, these shifts do not affect the main conclusions derived from Fig. 3c: in all sensitivity analyses, growth rate changes remain the main driver of global forest biomass C sink with total net C-sinks in countries where increasing growth rate is the main driver (including both solid and hatched countries) ranging from 12.1 to 21.1 GtC, while afforestation always holds the second place of global C sink driver (total net C-sinks in countries where afforestation is the main driver ranging from 2.4 to 7.7 GtC). Similarly, total net C sources by countries where deforestation is the main driver range from −21.9 to −14.0 GtC, thus highlighting that deforestation would by far remain the main driver of forest biomass C emissions across all sensitivity analyses.Implications for forest-based solutionsOur results allow to identify major mechanisms behind observed forest biomass C changes that are immediately relevant for forest-based climate-change-mitigation strategies. We show that deforestation, increasing harvest, and burnt area have driven the net C emissions from forest biomass over the last three decades. Deforestation is the dominant driver, corroborating that protection from deforestation is indispensable1,11,38. On the other hand, forest growth rate is identified as the major driver counteracting net C emissions (Fig. 2a, d). In fact, most of the temperate and boreal countries, with the noteworthy exception of Canada, fall under a type in which enhanced forest growth rate is the major driver of a net C sink (Fig. 3b). Besides, even countries dominated by deforestation in the tropics show significant increases in growth rate (Figs. 2d and 4). These results highlight that enhanced growth rate, rather than reforestation, is the main driver counteracting biomass C emissions in 1990–2020.Fig. 4: Change in forest growth rate and its effects on global carbon stocks.The diagrams show national forest growth rate changes (y-axis) scaled along the cumulated size of the carbon stock in 1990 (x-axis). The area between the graph and the x-axis indicates the C-stock change due to growth rate for total (a), primary (b), and managed forests (c) (see Supplementary Figs. 8–10 for results from sensitivity analyses).Full size imageThese increases in forest growth rate may arise from diverse processes, including climatic and land-use drivers. On the one hand, several studies highlight the effects of environmental drivers—such as warming, atmospheric carbon dioxide (CO2), and nitrogen (N) fertilization1,6,8,11,21,39—on the terrestrial C sink. On the other hand, changes in forest growth rate can also be driven by shifts in forest management practices, such as tree species selection, forest recovery from past degradation and lesser litter grazing12,40,41. Advancing the understanding of the underlying processes of forest growth rate change is key for forging climate-change-mitigation strategies, but it is not straightforward to isolate climatic (e.g., altered CO2 concentration or temperature) from land-use drivers (e.g., non-timber forest uses such as grazing)42. Still, a comparison of trajectories in primary and managed forest growth rate change based on our results allows to derive insights into the interplay of these different drivers (Fig. 4 and Supplementary Fig. 3). From the fact that only 11% of primary forest carbon stocks show declining growth rate trends (Fig. 4c) while a relatively larger carbon stock in managed forest (22%) is affected by declining growth rate trends (Fig. 4b), we can infer that in overall terms—and assuming primary and managed forests of a given country to be similarly affected by climatic drivers —land use is likely to exert a degrading effect on growth rate dynamics. Nevertheless, some countries reveal declining growth rate in primary forest but increasing growth rate in managed forest, thus suggesting that forest management may have an improving effect on forest growth rate in those countries (e.g., USA, Fig. 4b, c, see also Supplementary Fig. 4). In overall terms, this result suggests that globally a reduction of forest use may have the potential to enhance growth rate, thus corroborating previous findings by Quesada et al.14. However, these interpretations warrant a caveat that primary versus managed forest growth rate changes are derived from the FRA data and a state-of-the-art of the literature on changes in primary forest density (see “Methods” section and Supplementary Note 2), the latter being associated with higher uncertainties although the corresponding sensitivity analysis testing suggests these uncertainties to have little impact on the figures displayed here (see Supplementary Tables 2-3 and Supplementary Figs. 5–10).Independent of their origin (management or climate driven), the future trajectories of this driver, forest growth rate, is subject to large uncertainties43,44,45. Research suggests that increasing forest growth rate is a transient phenomenon and might be discontinued in the future46. For instance, several recent studies have pointed toward the saturating effect of CO2 fertilization, which is suspected to be a key process underlying vegetation greening and ensuing thickening21, the risk of increasing mortality and slower growth rate following increasing drought6,47,48, temperature49, and natural disturbances such as insect outbreaks50,51. Even more recently, Duffy et al.52 showed that, in the near-future, temperature increases from business-as-usual trajectories of climate change shall result in a severe reduction, and possibly a reversal, of the terrestrial C sink, despite the remaining unknowns.Therefore, we conclude that, while increasing forest growth rate is the dominant driver counteracting the global net C emissions from forest biomass in the past three decades, it is against a precautionary principle to forge climate strategies that rely on a continuous net C sink effect from the same processes in the future. By contrast, our results suggest that reducing wood harvest (Fig. 2g) and halting deforestation (Fig. 2c) are key strategies to address the challenge of climate-change mitigation. In this context, increasing forest harvest volumes—a strategy often promoted in the course of climate-change-mitigation efforts embraced as the “bioeconomy”—appears to have critical unintended side-effects, despite the potential of wood for substituting some emissions-intensive products and processes53,54,55: by not only reducing the carbon sink function in forests, but also accelerating the overall C turnover rates through rejuvenation of forests and transfer to harvested wood products of lifetimes shorter than those of old-growth forests56,57,58, such strategies result in a critical loss of C sink capacity. Overall, our results plead for a double strategy to enable future forest-based solutions for climate-change mitigation: in the Global South, ending deforestation is the main priority to reverse the net C source toward a net C sink, while in the Global North, lowering wood harvest has the strongest potential to immediately enhance the C sink in forest biomass. More

1.BirdLife International. State of the World’s Birds: Taking the Pulse of the Planet (BirdLife International, 2018).
Google Scholar 
2.Dias, M. P. et al. Threats to seabirds: A global assessment. Biol. Conserv. 237, 525–537 (2019).Article 

Google Scholar 
3.Gales, R. in Albatross Biology and Conservation (eds Robertson, G. & Gales, R.) 20–45 (Surrey Beatty and Sons, Chipping Norton, 1998).4.Gales, R., Brothers, N. & Reid, T. Seabird mortality in the Japanese tuna longline fishery around Australia, 1988–1995. Biol. Conserv. 86, 37–56 (1998).Article 

Google Scholar 
5.Anderson, O. R. et al. Global seabird bycatch in longline fisheries. Endanger. Species Res. 14, 91–106 (2011).Article 

Google Scholar 
6.Warham, J. The Petrels: Their Ecology and Breeding Systems (Academic Press, 1990).
Google Scholar 
7.Warham, J. The Behaviour, Population Biology and Physiology of the Petrels (Academic Press, 1996).
Google Scholar 
8.Dietrich, K. S., Parrish, J. K. & Melvin, E. F. Understanding and addressing seabird bycatch in Alaska demersal longline fisheries. Biol. Conserv. 142, 2642–2656 (2009).Article 

Google Scholar 
9.Zhou, C., Jiao, Y. & Browder, J. A. Seabird bycatch vulnerability to pelagic longline fisheries: Ecological traits matter. Aquat. Conserv.: Mar. Freshw. Ecosyst. https://doi.org/10.1002/aqc.3066 (2019).Article 

Google Scholar 
10.Brothers, N. The incidental catch of seabirds by longline fisheries: Worldwide review and technical guidelines for mitigation. FAO Fish. Circ. 937, 1–100 (1999).ADS 

Google Scholar 
11.Gilman, E. Integrated management to address the incidental mortality of seabirds in longline fisheries. Aquat. Conserv.: Mar. Freshw. Ecosyst. 11, 391–414 (2001).Article 

Google Scholar 
12.Li, Y. & Jiao, Y. Modeling spatial patterns of rare species using eigenfunction-based spatial filters: An example of modified delta model for zero-inflated data. Ecol. Model. 299, 51–63 (2015).Article 

Google Scholar 
13.Li, Y., Jiao, Y. & Browder, J. A. Assessment of seabird bycatch in the US Atlantic pelagic longline fishery, with an extra exploration on modeling spatial variation. ICES J. Mar. Sci. 73, 2687–2694 (2016).Article 

Google Scholar 
14.Brothers, N. Albatross mortality and associated bait loss in the Japanese longline fishery in the Southern Ocean. Biol. Conserv. 55, 255–268 (1991).Article 

Google Scholar 
15.Croxall, J. P. & Nicol, S. Management of Southern Ocean fisheries: Global forces and future sustainability. Antarct. Sci. 16, 569–584 (2004).ADS 
Article 

Google Scholar 
16.Løkkeborg, S. Best practices to mitigate seabird bycatch in longline, trawl and gillnet fisheries – efficiency and practical applicability. Mar. Ecol. Prog. Ser. 435, 285–303 (2011).ADS 
Article 

Google Scholar 
17.Beerkircher, L. R., Brown, C. J., Abercrombie, D. L. & Lee, D. W. Overview of the SEFSC pelagic observer program in the Northwest Atlantic from 1992–2002. ICCAT CVSP 58, 1729–1748 (2005).
Google Scholar 
18.Diaz, G. A., Beerkircher, L. R. & Restrepo, V. R. Description of the US pelagic observer program (POP). ICCAT CVSP 64, 2415–2426 (2009).
Google Scholar 
19.Lee, D. W. & Brown, C. J. SEFSC pelagic observer program data summary for 1992–1996. US Department of Commerce, National Oceanic and Atmospheric Administration, National Marine Fisheries Service, Southeast Fisheries Science Center (1998).20.Lo, N. C., Jacobson, L. D. & Squire, J. L. Indices of relative abundance from fish spotter data based on delta-lognornial models. Can. J. Fish. Aquat. Sci. 49, 2515–2526 (1992).Article 

Google Scholar 
21.Martin, T. G. et al. Zero tolerance ecology: improving ecological inference by modelling the source of zero observations. Ecol. Lett. 8, 1235–1246 (2005).PubMed 
Article 
PubMed Central 

Google Scholar 
22.Winter, A., Jiao, Y. & Browder, J. A. Modeling low rates of seabird bycatch in the US Atlantic longline fishery. Waterbirds 34, 289–303 (2011).Article 

Google Scholar 
23.Cortés, V., Arcos, J. M. & González-Solís, J. Seabirds and demersal longliners in the northwestern Mediterranean: Factors driving their interactions and bycatch rates. Mar. Ecol. Prog. Ser. 565, 1–16 (2017).ADS 
Article 

Google Scholar 
24.Bi, R., Jiao, Y., Zhou, C. & Hallerman, E. M. A Bayesian spatiotemporal approach to inform management unit appropriateness. Can. J. Fish. Aquat. Sci. 76, 217–237 (2018).Article 

Google Scholar 
25.Wikle, C. K. Hierarchical models in environmental science. Int. Stat. Rev. 71, 181–199 (2003).MATH 
Article 

Google Scholar 
26.Cressie, N., Calder, C. A., Clark, J. S., Hoef, J. M. V. & Wikle, C. K. Accounting for uncertainty in ecological analysis: The strengths and limitations of hierarchical statistical modeling. Ecol. Appl. 19, 553–570 (2009).PubMed 
Article 
PubMed Central 

Google Scholar 
27.Banerjee, S., Carlin, B. P. & Gelfand, A. E. Hierarchical Modeling and Analysis for Spatial Data (CRC Press, 2014).MATH 
Book 

Google Scholar 
28.Besag, J., York, J. & Mollié, A. Bayesian image restoration, with two applications in spatial statistics. Ann. Inst. Stat. Math. 43, 1–20 (1991).MathSciNet 
MATH 
Article 

Google Scholar 
29.Rue, H. & Held, L. Gaussian Markov Random Fields: Theory and Applications (CRC Press, 2005).MATH 
Book 

Google Scholar 
30.Rue, H., Martino, S. & Chopin, N. Approximate Bayesian inference for latent Gaussian models by using integrated nested Laplace approximations. J. R. Stat. Soc. Ser. B Stat. Methodol. 71, 319–392 (2009).MathSciNet 
MATH 
Article 

Google Scholar 
31.Held, L., Schrödle, B. & Rue, H. in Statistical Modelling and Regression Structures (eds Kneib, T. & Tutz, G.) 91–110 (Springer, Berlin Heidelberg, 2010).32.Lindgren, F., Rue, H. & Lindström, J. An explicit link between Gaussian fields and Gaussian Markov random fields: The stochastic partial differential equation approach. J. R. Stat. Soc. Ser. B Stat. Methodol. 73, 423–498 (2011).MathSciNet 
MATH 
Article 

Google Scholar 
33.Bakka, H., Vanhatalo, J., Illian, J. B., Simpson, D. & Rue, H. Non-stationary Gaussian models with physical barriers. Spat. Stat. 29, 268–288 (2019).MathSciNet 
Article 

Google Scholar 
34.NCAR (National Center for Atmospheric Research). The Climate Data Guide: Hurrell North Atlantic Oscillation (NAO) Index (PC-based). Available from: https://climatedataguide.ucar.edu/climate-data/hurrell-north-atlantic-oscillation-nao-index-pc-based. Retrieved: March 1, 2019.35.ESRL (Earth Science Research Laboratory, NOAA). Climate timeseries: AMO (Atlantic Multidecadal Oscillation) Index. Available from: http://www.esrl.noaa.gov/psd/data/timeseries/AMO/. Retrieved: March 1, 2019.36.Spiegelhalter, D. J., Best, N. G., Carlin, B. P. & Van Der Linde, A. Bayesian measures of model complexity and fit. J. R. Stat. Soc. Ser. B Stat. Methodol. 64, 583–639 (2002).MathSciNet 
MATH 
Article 

Google Scholar 
37.Watanabe, S. Asymptotic equivalence of Bayes cross validation and widely applicable information criterion in singular learning theory. J. Mach. Learn. Res. 11, 3571–3594 (2010).MathSciNet 
MATH 

Google Scholar 
38.Shumway, R. H. & Stoffer, D. S. Time Series Analysis and Its Applications (Springer, 2011).MATH 
Book 

Google Scholar 
39.Bi, R., Jiao, Y., Bakka, H. & Browder, J. A. Long-term climate ocean oscillations inform seabird bycatch from pelagic longline fishery. ICES J. Mar. Sci. 77, 668–679 (2020).Article 

Google Scholar 
40.Lear, W. H. History of fisheries in the Northwest Atlantic: The 500 year perspective. J. Northwest Atl. Fish. Sci. 23, 41–73 (1998).Article 

Google Scholar 
41.Veit, R. R., Goyert, H. F., White, T. P., Martin, M. C., Manne, L. L. & Gilbert, A. Pelagic Seabirds off the East Coast of the United States 2008–2013. US Dept. of the Interior, Bureau of Ocean Energy Management, Office of Renewable Energy Programs, Sterling, VA. OCS Study BOEM, 24, 186 (2015).42.Harrison, P. Seabirds, an identification guide (Houghton Mifflin, 1983).
Google Scholar 
43.Onley, D. & Scofield, P. Albatrosses, petrels and shearwaters of the world (Princeton University Press, 2013).
Google Scholar 
44.Gladics, A. J. et al. Fishery-specific solutions to seabird bycatch in the U.S. West Coast sablefish fishery. Fish. Res. 196, 85–95 (2017).Article 

Google Scholar 
45.Grieve, B. D., Hare, J. A. & Saba, V. S. Projecting the effects of climate change on Calanus finmarchicus distribution within the U.S. Northeast Continental Shelf. Sci. Rep. 7, 6264 (2017).ADS 
PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
46.Petersen, S. L., Honig, M. B., Ryan, P. G. & Underhill, L. G. Seabird bycatch in the pelagic longline fishery off southern Africa. Afr. J. Mar. Sci. 31, 191–204 (2009).Article 

Google Scholar 
47.Arcos, J. M. & Oro, D. Significance of fisheries discards for a threatened Mediterranean seabird, the Balearic shearwater Puffinus mauretanicus. Mar. Ecol. Prog. Ser. 239, 209–220 (2002).ADS 
Article 

Google Scholar 
48.Furness, R., Edwards, A. & Oro, D. Influence of management practices and of scavenging seabirds on availability of fisheries discards to benthic scavengers. Mar. Ecol. Prog. Ser. 350, 235–244 (2007).ADS 
Article 

Google Scholar 
49.Grémillet, D. et al. A junk-food hypothesis for gannets feeding on fishery waste. Proc. Biol. Sci. 275, 1149–1156 (2008).PubMed 
PubMed Central 

Google Scholar 
50.Skov, H. & Durinck, J. Seabird attraction to fishing vessels is a local process. Mar. Ecol. Prog. Ser. 214, 289–298 (2001).ADS 
Article 

Google Scholar 
51.Chapman, D. C., Barth, J. A., Beardsley, R. C. & Fairbanks, R. G. On the continuity of mean flow between the Scotian Shelf and the Middle Atlantic Bight. J. Phys. Oceanogr. 16, 758–772 (1986).ADS 
Article 

Google Scholar 
52.Steimle, F. W. & Zetlin, C. Reef habitats in the middle Atlantic bight: Abundance, distribution, associated biological communities, and fishery resource use. Mar. Fish. Rev. 62, 24–42 (2000).
Google Scholar 
53.Lee, D. S. Pelagic seabirds and the proposed exploration for fossil fuels off North Carolina: A test for conservation efforts of a vulnerable international resource. J. Elisha Mitchell Sci. Soc. 115, 294–315 (1999).
Google Scholar 
54.Kai, E. T. et al. Top marine predators track Lagrangian coherent structures. Proc. Natl. Acad. Sci. 106, 8245–8250 (2009).ADS 
CAS 
Article 

Google Scholar 
55.Li, Y., Browder, J. A. & Jiao, Y. Hook effects on seabird bycatch in the United States Atlantic pelagic longline fishery. Bull. Mar. Sci. 88, 559–569 (2012).Article 

Google Scholar 
56.Taylor, A. H. & Stephens, J. A. The North Atlantic Oscillation and the latitude of the Gulf Stream. Tellus 50, 134–142 (1998).Article 

Google Scholar 
57.Hobday, A. J., Hartog, J. R., Spillman, C. M. & Alves, O. Seasonal forecasting of tuna habitat for dynamic spatial management. Can. J. Fish. Aquat. Sci. 68, 898–911 (2011).Article 

Google Scholar 
58.FAO (Food and Agriculture Organization of the United Nations). Guidelines to reduce sea turtle mortality in fishing operations. FAO Technical Guidelines for Responsible Fisheries Prepared by Gilman, E., Bianchi, G. FAO: Rome. ISBN 978-92-106226-5 (2009).59.Bethoney, N. D., Schondelmeier, B. P., Kneebone, J. & Hoffman, W. S. Bridges to best management: Effects of a voluntary bycatch avoidance program in a mid-water trawl fishery. Mar. Policy 83, 172–178 (2017).Article 

Google Scholar 
60.Lindgren, F. & Rue, H. Bayesian spatial modelling with R-INLA. J. Stat. Softw. 63, 19 (2015).Article 

Google Scholar 
61.R Core Team. R: A language and environment for statistical computing. R Foundation for Statistical Computing, Vienna, Austria. https://www.R-project.org/ (2019).62.Rue, H., Martino, S. & Lindgren, F. INLA: Functions which allow to perform a full Bayesian analysis of structured (geo-)additive models using integrated nested Laplace approximation. R package version 0.0., GNU General Public License, version 3 (2009).63.Simpson, D., Rue, H., Riebler, A., Martins, T. G. & Sørbye, S. H. Penalising model component complexity: A principled, practical approach to constructing priors (with discussion). Stat. Sci. 32, 1–28 (2017).MATH 

Google Scholar 
64.Fuglstad, G.-A., Simpson, D., Lindgren, F. & Rue, H. Constructing priors that penalize the complexity of Gaussian random fields. J. Am. Stat. Assoc. 114, 445–452 (2018).MathSciNet 
MATH 
Article 
CAS 

Google Scholar 
65.Plummer, M. Penalized loss functions for Bayesian model comparison. Biostat. 9, 523–539 (2008).MATH 
Article 

Google Scholar 
66.Gelman, A., Hwang, J. & Vehtari, A. Understanding predictive information criteria for Bayesian models. Stat. Comput. 24, 997–1016 (2014).MathSciNet 
MATH 
Article 

Google Scholar  More

1.Krebs, C. J. Two paradigms of population regulation. Wildl. Res. 22, 1–10 (1995).MathSciNet 
Article 

Google Scholar 
2.McLaren, I. A. Natural Regulation of Animal Populations (Routledge, 2017).Book 

Google Scholar 
3.Deffner, D. & McElreath, R. The importance of life history and population regulation for the evolution of social learning. Philos. Trans. R. Soc. B 375, 20190492 (2020).Article 

Google Scholar 
4.Leão, S. M., Pianka, E. R. & Pelegrin, N. Is there evidence for population regulation in amphibians and reptiles? J. Herpetol. 52, 28–33 (2018).Article 

Google Scholar 
5.Hanski, I. A. Density dependence, regulation and variability in animal populations. Philos. Trans. R. Soc. B 330, 141–150 (1990).ADS 
Article 

Google Scholar 
6.Reznick, D., Bryant, M. J. & Bashey, F. r-and K-selection revisited: The role of population regulation in life-history evolution. Ecology 83, 1509–1520 (2002).Article 

Google Scholar 
7.Stenseth, N. C., Falck, W., Bjørnstad, O. N. & Krebs, C. J. Population regulation in snowshoe hare and Canadian lynx: Asymmetric food web configurations between hare and lynx. Proc. Natl Acad. Sci. USA 94, 5147–5152 (1997).ADS 
CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
8.Lande, R. et al. Estimating density dependence from population time series using demographic theory and life-history data. Am. Nat. 159, 321–337 (2002).CAS 
PubMed 
Article 

Google Scholar 
9.Ferguson, G. W., Bohlen, C. H. & Woolley, H. P. Sceloporus undulatus: Comparative life history and regulation of a Kansas population. Ecology 61, 313–322 (1980).Article 

Google Scholar 
10.Andrews, R. M. Population stability of a tropical lizard. Ecology 72, 1204–1217 (1991).Article 

Google Scholar 
11.Tinkle, D. W., Dunham, A. E. & Congdon, J. D. Life history and demographic variation in the lizard Sceloporus graciosus: A long-term study. Ecology 74, 2413–2429 (1993).Article 

Google Scholar 
12.Altwegg, R., Dummermuth, S., Anholt, B. R. & Flatt, T. Winter weather affects asp viper Vipera aspis population dynamics through susceptible juveniles. Oikos 110, 55–66 (2005).Article 

Google Scholar 
13.Madsen, T. & Shine, R. Rain, fish and snakes: Climatically driven population dynamics of Arafura filesnakes in tropical Australia. Oecologia 124, 208–215 (2000).ADS 
CAS 
PubMed 
Article 

Google Scholar 
14.Madsen, T., Ujvari, B., Shine, R. & Olsson, M. Rain, rats and pythons: Climate-driven population dynamics of predators and prey in tropical Australia. Austral Ecol. 31, 30–37 (2006).Article 

Google Scholar 
15.Brown, G. P. & Shine, R. Rain, prey and predators: Climatically driven shifts in frog abundance modify reproductive allometry in a tropical snake. Oecologia 154, 361–368 (2007).ADS 
PubMed 
Article 

Google Scholar 
16.Brown, G. P., Ujvari, B., Madsen, T. & Shine, R. Invader impact clarifies the roles of top-down and bottom-up effects on tropical snake populations. Funct. Ecol. 27, 351–361 (2013).Article 

Google Scholar 
17.Massot, M., Clobert, J., Pilorge, T., Lecomte, J. & Barbault, R. Density dependence in the common lizard: Demographic consequences of a density manipulation. Ecology 73, 1742–1756 (1992).Article 

Google Scholar 
18.Fordham, D. A., Georges, A. & Brook, B. W. Experimental evidence for density-dependent responses to mortality of snake-necked turtles. Oecologia 159, 271–281 (2009).ADS 
PubMed 
Article 

Google Scholar 
19.Burns, G. & Heatwole, H. Home range and habitat use of the olive sea snake, Aipysurus laevis, on the Great Barrier Reef, Australia. J. Herpetol. 32, 350–358 (1998).Article 

Google Scholar 
20.Ward, T. M. Age structures and reproductive patterns of two species of sea snake, Lapemis hardwickii Grey (1836) and Hydrophis elegans (Grey 1842), incidentally captured by prawn trawlers in northern Australia. Mar. Freshw. Res. 52, 193–203 (2001).Article 

Google Scholar 
21.Dennis, B. & Ponciano, J. M. Density-dependent state-space model for population-abundance data with unequal time intervals. Ecology 95, 2069–2076 (2014).PubMed 
Article 

Google Scholar 
22.Bonnet, X., Naulleau, G. & Shine, R. The dangers of leaving home: Dispersal and mortality in snakes. Biol. Conserv. 89, 39–50 (1999).Article 

Google Scholar 
23.Yacelga, M., Cayot, L. J. & Jaramillo, A. Dispersal of neonatal Galápagos marine iguanas Amblyrhynchus cristatus from their nesting zone: Natural history and conservation implications. Herpetol. Conserv. Biol. 7, 470–480 (2012).
Google Scholar 
24.Lukoschek, V. & Shine, R. Sea snakes rarely venture far from home. Ecol. Evol. 2, 1113–1121 (2012).PubMed 
PubMed Central 
Article 

Google Scholar 
25.Shine, R., Cogger, H. G., Reed, R. R., Shetty, S. & Bonnet, X. Aquatic and terrestrial locomotor speeds of amphibious sea-snakes (Serpentes, Laticaudidae). J. Zool. 259, 261–268 (2003).Article 

Google Scholar 
26.Forsman, A. Body size and net energy gain in gape-limited predators: A model. J. Herpetol. 30, 307–319 (1996).Article 

Google Scholar 
27.Shine, R., LeMaster, M. P., Moore, I. T., Olsson, M. M. & Mason, R. T. Bumpus in the snake den: Effects of sex, size and body condition on mortality in red-sided garter snakes. Evolution 55, 598–604 (2001).CAS 
PubMed 
Article 

Google Scholar 
28.Sherratt, E., Rasmussen, A. R. & Sanders, K. L. Trophic specialization drives morphological evolution in sea snakes. R. Soc. Open Sci. 5, 172141 (2018).ADS 
PubMed 
PubMed Central 
Article 

Google Scholar 
29.Lemen, C. A. & Voris, H. K. A comparison of reproductive strategies among marine snakes. J. Anim. Ecol. 50, 89–101 (1981).Article 

Google Scholar 
30.Shine, R., Shine, T. & Shine, B. Intraspecific habitat partitioning by the sea snake Emydocephalus annulatus (Serpentes, Hydrophiidae): The effects of sex, body size, and colour pattern. Biol. J. Linn. Soc. 80, 1–10 (2003).Article 

Google Scholar 
31.Heatwole, H., Grech, A., Monahan, J. F., King, S. & Marsh, H. Thermal biology of sea snakes and sea kraits. Integr. Comp. Biol. 52, 257–273 (2012).PubMed 
Article 

Google Scholar 
32.Van Dyke, J. U., Beaupre, S. J. & Kreider, D. L. Snakes allocate amino acids acquired during vitellogenesis to offspring: Are capital and income breeding consequences of variable foraging success? Biol. J. Linn. Soc. 106, 390–404 (2012).Article 

Google Scholar 
33.Masunaga, G., Matsuura, R., Yoshino, T. & Ota, H. Reproductive biology of the viviparous sea snake Emydocephalus ijimae (Reptilia: Elapidae: Hydrophiinae) under a seasonal environment in the Northern Hemisphere. Herpetol. J. 13, 113–119 (2003).
Google Scholar 
34.Goiran, C., Dubey, S. & Shine, R. Effects of season, sex and body size on the feeding ecology of turtle-headed sea snakes (Emydocephalus annulatus) on IndoPacific inshore coral reefs. Coral Reefs 32, 527–538 (2013).ADS 
Article 

Google Scholar 
35.Phillips, B. L. The evolution of growth rates on an expanding range edge. Biol. Lett. 5, 802–804 (2009).PubMed 
PubMed Central 
Article 

Google Scholar 
36.Shine, R., Shine, T. G., Brown, G. P. & Goiran, C. Life history traits of the sea snake Emydocephalus annulatus, based on a 17-yr study. Coral Reefs 39, 1407–1414 (2020).Article 

Google Scholar 
37.Goiran, C. & Shine, R. Decline in sea snake abundance on a protected coral reef system in the New Caledonian Lagoon. Coral Reefs 32, 281–284 (2013).ADS 
Article 

Google Scholar 
38.Somaweera, R. et al. Pinpointing drivers of extirpation in sea snakes: A synthesis of evidence from Ashmore Reef. Front. Mar. Sci. 8, 658756 (2021).Article 

Google Scholar 
39.Udyawer, V. et al. Future directions in the research and management of marine snakes. Front. Mar. Sci. 5, 399 (2018).Article 

Google Scholar 
40.Udyawer, V., Cappo, M., Simpfendorfer, C. A., Heupel, M. R. & Lukoschek, V. Distribution of sea snakes in the Great Barrier Reef Marine Park: Observations from 10 yrs of baited remote underwater video station (BRUVS) sampling. Coral Reefs 33, 777–791 (2014).ADS 
Article 

Google Scholar 
41.Udyawer, V., Goiran, C. & Shine, R. Peaceful coexistence between people and deadly wildlife: Why are recreational users of the ocean so rarely bitten by sea snakes? People Nat. 3, 335–346 (2021).Article 

Google Scholar 
42.Shine, R., Goiran, C., Shine, T., Fauvel, T. & Brischoux, F. Phenotypic divergence between seasnake (Emydocephalus annulatus) populations from adjacent bays of the New Caledonian Lagoon. Biol. J. Linn. Soc. 107, 824–832 (2012).Article 

Google Scholar 
43.Goiran, C., Brown, G. P. & Shine, R. Niche partitioning within a population of sea snakes is constrained by ambient thermal homogeneity and small prey size. Biol. J. Linn. Soc. 129, 644–651 (2020).Article 

Google Scholar 
44.Li, M., Fry, B. G. & Kini, R. M. Eggs-only diet: Its implications for the toxin profile changes and ecology of the marbled sea snake (Aipysurus eydouxii). J. Mol. Evol. 60, 81–89 (2005).ADS 
CAS 
PubMed 
Article 

Google Scholar 
45.Heatwole, H. Predation on sea snakes. In The Biology of Sea Snakes (ed. Dunson, W. A.) 233–250 (University Park Press, 1975).
Google Scholar 
46.Rancurel, P. & Intes, A. L. requin tigre, Galeocerdo cuvieri Lacepede, des eaux neocaledoniennes examen des contenus stomacaux. Tethys 10, 195–199 (1982).
Google Scholar 
47.Ineich, I. & Laboute, P. Les Serpents Marins de Nouvelle-Calédonie (IRD éditions, 2002).
Google Scholar 
48.Masunaga, G., Kosuge, T., Asai, N. & Ota, H. Shark predation of sea snakes (Reptilia: Elapidae) in the shallow waters around the Yaeyama Islands of the southern Ryukyus, Japan. Mar. Biodivers. Rec. 1, e96 (2008).Article 

Google Scholar 
49.Wirsing, A. J. & Heithaus, M. R. Olive-headed sea snakes Disteria major shift seagrass microhabitats to avoid shark predation. Mar. Ecol. Progr. Ser. 387, 287–293 (2009).ADS 
Article 

Google Scholar 
50.White, G. C. & Burnham, K. P. Program MARK: Survival estimation from populations of marked animals. Bird Study 46, S120–S139 (1999).Article 

Google Scholar 
51.Brischoux, F., Rolland, V., Bonnet, X., Caillaud, M. & Shine, R. Effects of oceanic salinity on body condition in sea snakes. Integr. Comp. Biol. 52, 235–244 (2012).PubMed 
Article 

Google Scholar 
52.Lovich, J. E. & Gibbons, J. W. A review of techniques for quantifying sexual size dimorphism. Growth Dev. Aging 56, 269–269 (1992).CAS 
PubMed 

Google Scholar 
53.Dennis, B., Ponciano, J. M., Lele, S. R., Taper, M. L. & Staples, D. F. Estimating density dependence, process noise, and observation error. Ecol. Monogr. 76, 323–341 (2006).Article 

Google Scholar  More