Ecology

Phenological and climatic data
We used data from the Pan European Phenology Project (PEP725)29, an open-access database with long-term plant phenological observations across 25 European countries (http://www.pep725.eu/). The regional/national network partners of PEP725 are following a consistent guideline for phenological observations30 and prepare the data for submission to the PEP725 database curators29. We selected 30 species for which sufficient observational data were available: 21 deciduous broadleaved trees or shrubs, 6 herbaceous perennials, 2 evergreen coniferous trees, and 1 deciduous coniferous tree (Supplementary Table 3). Particularly, our data set included one fruit tree (Prunus avium) and one nut tree (Corylus avellana) since some of the chilling models are specifically developed for fruit and nut trees. A total of 2,493,644 individual records from 15,533 phenological stations were used. The stations were mainly distributed in moderate climates in Central Europe (Supplementary Fig. 3). Four spring events based on the BBCH code were investigated: BBCH 10, 11, 60, and 69, representing first leaves separated, first leaves unfolded, first flowers open, and end of flowering, respectively31.
We used the E-OBS v19.0eHOM data set32 with a spatial resolution of 0.1 × 0.1° for 1950–2018 for calculating CA and HR of the in situ phenological records. This data set is provided by the European Climate Assessment & Data set project and includes homogenized series of daily mean, minimum, and maximum temperatures. We also use the daily maximum and minimum temperature data from the GHCN data set33 to assess the scale effect. The GHCN data set contains station-based measurements from over 90,000 land-based stations worldwide, but only parts of PEP725 stations match with the GHCN stations.
For future climatic data (2019–2099), we used daily minimum and maximum temperatures simulated by the HADGEM2-ES model (with a spatial resolution of 0.5 × 0.5°) under two climatic scenarios (RCP 4.5 and RCP 8.5). These data have been bias-corrected by applying the method used in the Inter-Sectoral Impact Model Intercomparison Project (ISIMIP)34, which were available on the ISIMIP server (https://esg.pik-potsdam.de/projects/isimip2b/).
Chilling models
We used 12 chilling models to measure the amount of chilling. One type of chilling model is based on several specific temperature thresholds. The most commonly used model, developed in the 1930s and 1940s for peach35, calculates the number of hours or days with temperatures 0 °C for calculating CA11,38. Models C1–C6 were developed based on various combinations of the upper and lower temperature limits. The rate of chilling was 1 for daily temperatures 5,{mathrm{or}},{T} < - 10} hfill end{array}} right.,$$ (2) $${mathop{rm{CU}}nolimits} _3 = left{ {begin{array}{*{20}{l}} 1 hfill & {0 le {mathop{T}nolimits} le 5} hfill \ 0 hfill & {{T} > 5,{mathrm{or}},{T} < 0} hfill end{array}} right.,$$ (3) $${mathrm{CU}}_{mathrm{4}} = left{ {begin{array}{*{20}{c}} 1 & {{mathop{T}nolimits} le 7} \ 0 & {{T} > 7} end{array}} right.,$$
(4)

$${mathop{rm{CU}}nolimits} _5 = left{ {begin{array}{*{20}{l}} 1 hfill & { – 10 le {mathop{T}nolimits} le 7} hfill \ 0 hfill & {{T} > 7,{mathrm{or}},{T} < - 10} hfill end{array}} right.,$$ (5) $${mathop{rm{CU}}nolimits} _6 = left{ {begin{array}{*{20}{l}} 1 hfill & {0 le {mathop{T}nolimits} le 7} hfill \ 0 hfill & {{T} > 7,{mathrm{or}},{T} < 0} hfill end{array}} right.,$$ (6) where CUi is the rate of chilling for Model Ci, and T is the daily mean temperature (°C). Model C7 is also known as the Utah Model39, which assigned different weights to different ranges of temperatures and was first used to measure the chilling requirements of peach (Eq. (7)). The Utah Model was modified to produce Model C8 (Eq. (8)), which removed the negative contributions of warm temperatures to accumulated chilling20. $${mathop{rm{CU}}nolimits} _7 = left{ {begin{array}{*{20}{l}} 0 hfill & {{mathop{T}nolimits} le 1.4} hfill \ {0.5} hfill & {1.4 < {mathop{T}nolimits} le 2.4} hfill \ 1 hfill & {2.4 < {mathop{T}nolimits} le 9.1} hfill \ {0.5} hfill & {9.1 < {mathop{T}nolimits} le 12.4} hfill \ 0 hfill & {12.4 < {mathop{T}nolimits} le 15.9} hfill \ { - 0.5} hfill & {15.9 < {mathop{T}nolimits} le 18} hfill \ { - 1} hfill & {{mathop{T}nolimits} > 18} hfill end{array}} right.,$$
(7)

$${mathop{rm{CU}}nolimits} _8 = left{ {begin{array}{*{20}{l}} 0 hfill & {{mathop{T}nolimits} le 1.4} hfill \ {0.5} hfill & {1.4 < {mathop{T}nolimits} le 2.4} hfill \ 1 hfill & {2.4 < {mathop{T}nolimits} le 9.1} hfill \ {0.5} hfill & {9.1 < {mathop{T}nolimits} le 12.4} hfill \ 0 hfill & {{mathop{T}nolimits} > 12.4} hfill end{array}} right.,$$
(8)

where CUi is the rate of chilling for Model Ci, and T is the daily mean temperature (°C).
Model C9 is a dynamic model developed for peach in Israel and South Africa40,41 and now adopted for apricot cultivars42. The most important characteristic of Model C9 was that a previous intermediate product affected the rate of chilling in the current hour or day. We did not provide equations for Model C9 for simplicity (see the equation in Luedeling et al.20).
Harrington et al.43 summarized published results for chilling units and constructed a chilling function based on a three-parameter Weibull distribution, coded as Model C10 (Eq. (9)). Model C11 has a triangular form, which was fitted by Hänninen44 using previous experimental results for Finnish birch seedlings (Eq. (10)). Zhang et al.45 recently fitted observational data to the triangular model for 24 plant species and found that a mean optimal chilling temperature of 0.2 °C and an upper limit of the chilling temperature of 6.9 °C were most effective. Model C12, therefore, uses the triangular form with parameters of 0.2 and 6.9 °C (Eq. (11)).

$${mathop{rm{CU}}nolimits} _{10} = left{ {begin{array}{*{20}{l}} 1 hfill & {2.5 < {mathop{T}nolimits} < 7.4} hfill \ 0 hfill & {{mathop{T}nolimits} < - 4.7,{mathrm{or}},{mathop{T}nolimits} > 16} hfill \ {3.13left( {frac{{{mathop{T}nolimits} + 4.66}}{{10.93}}} right)^{2.10}{mathop{rm{e}}nolimits} ^{ – left( {frac{{{mathop{T}nolimits} + 4.66}}{{10.93}}} right)^{3.10}}} hfill & {{mathrm{else}}} hfill end{array}} right.,$$
(9)

$${mathrm{CU}}_{11} = left{ {begin{array}{*{20}{l}} 0 hfill & {{mathop{T}nolimits} le – 3.4,{mathop{rm{or}}nolimits} ,{mathop{T}nolimits} ge 10.4} hfill \ {frac{{{mathop{T}nolimits} + 3.4}}{{5 + 3.4}}} hfill & { – 3.4 < {mathop{rm{T}}nolimits} le 5} hfill \ {frac{{{mathop{T}nolimits} - 10.4}}{{5 - 10.4}}} hfill & {5 < {mathop{T}nolimits} < 10.4} hfill end{array}} right.,$$ (10) $${mathop{rm{CU}}nolimits} _{12} = left{ {begin{array}{*{20}{l}} 0 hfill & {{mathop{T}nolimits} le - 6.5,{mathop{rm{or}}nolimits} ,{mathop{T}nolimits} ge 6.9} hfill \ {frac{{{T + }6.5}}{{6.9 - 0.2}}} hfill & { - 6.5 < {mathop{T}nolimits} le 0.2} hfill \ {frac{{6.9 - {T}}}{{6.9 - 0.2}}} hfill & {0.2 < {T} < 6.9} hfill end{array}} right.,$$ (11) where CUi is the rate of chilling for Model Ci, and T is the daily mean temperature in °C. Forcing models Forcing models were used to measure HR for the spring events of plants. The GDD model is the most commonly used forcing model, which assumes that the rate of forcing is linearly correlated with temperature if the temperature is above a particular threshold. We mainly used Model F1 (Eq. (12)), which adopts a temperature threshold of 0 °C46,47,48, for examining the relationship between CA and HR: $${mathrm{FU}}_{mathrm{1}} = {mathrm{max}}({T},0),$$ (12) where FU1 is the rate of forcing for Model F1, and T is the daily mean temperature (°C). We also validated the chilling models by correlating them with HR based on seven other forcing models to assess the impact of the choice of forcing model on the results. Model F2 (Eq. (13)) is also a GDD model but has a temperature threshold of 5 °C12,18: $${mathrm{FU}}_{mathrm{2}} = {mathrm{max}}(T - 5,0),$$ (13) where FU2 is the rate of forcing for Model F2, and T is the daily mean temperature (°C). Piao et al.48 found that leaf onset in the Northern Hemisphere was triggered more by daytime than nighttime temperature. They thus proposed a GDD model using maximum instead of mean temperature. Models F3 (Eq. (14)) and F4 (Eq. (15)) are based on maximum temperature with thresholds of 0 and 5 °C, respectively. $${mathop{rm{FU}}nolimits} _3 = max ({mathop{T}nolimits} _{rm{max}},0),$$ (14) $${mathop{rm{FU}}nolimits} _4 = max ({mathop{T}nolimits} _{rm{max}} - 5,0),$$ (15) where FUi is the rate of forcing for Model Fi, and Tmax is the daily maximum temperature (°C). A recent experiment demonstrated that the impact of daytime temperature on leaf unfolding for temperate trees was approximately threefold higher than the impact of nighttime temperature49. Model F5 (Eq. 16) thus uses two parameters (0.75 and 0.25) to weigh the impact of daytime and nighttime temperatures on HR. $${mathop{rm{FU}}nolimits} _5 = 0.75 times max ({mathop{T}nolimits} _{rm{max}} - 5,0){mathrm{ + }}0.25 times max ({mathop{T}nolimits} _{rm{min}} - 5,0),$$ (16) where FU5 is the rate of forcing for Model F5. Tmax and Tmin are the daily maximum and minimum temperatures (°C), respectively. Many studies have suggested that the rate of forcing followed a logistic function of temperature44,50. Model F6 (Eq. (17)) uses a logistic function proposed by Hänninen44, and Model F7 (Eq. (18)) uses another logistic function proposed by Harrington et al.43. $${mathop{rm{FU}}nolimits} _6 = left{ {begin{array}{*{20}{l}} {frac{{28.4}}{{1 + {mathop{rm{e}}nolimits} ^{ - 0.185({mathop{T}nolimits} - 18.5)}}}} hfill & {{mathop{T}nolimits} > 0} hfill \ 0 hfill & {{mathop{rm{else}}nolimits} } hfill end{array}} right.,$$
(17)

$${mathop{rm{FU}}nolimits} _7 = frac{1}{{1 + {mathop{rm{e}}nolimits} ^{ – 0.47{mathop{T}nolimits} {mathrm{ + 6}}{mathrm{.49}}}}},$$
(18)

where FUi is the rate of forcing for Model Fi, and T is the daily mean temperature (°C).
Model F8 is a growing degree hour (GDH) model, where species have an optimum temperature for growth and where temperatures above or below that optimum have a smaller impact51. Model F8 (Eq. (19)) was first designed for calculating HR at hourly intervals, but we applied it at a daily interval. The stress factor in the original GDH model was ignored, because we assumed that the plants were not under other stresses.

$${mathrm{FU}}_{mathrm{8}}{mathrm{ = }}left{ {begin{array}{*{20}{l}} {mathrm{0}} hfill & {{T < }T_{mathrm{L}},{mathrm{or}},{T > }T_{mathrm{c}}} hfill \ {frac{{T_{mathrm{u}} – T_{mathrm{L}}}}{{mathrm{2}}}left( {{mathrm{1 + cos}}left( {pi {mathrm{ + }}pi frac{{{T} – T_{mathrm{L}}}}{{T_{mathrm{u}} – T_{mathrm{L}}}}} right)} right)} hfill & {T_{mathrm{L}} ge {T} ge T_{mathrm{u}}} hfill \ {left( {T_{mathrm{u}} – T_{mathrm{L}}} right)left( {{mathrm{1 + cos}}left( {frac{{uppi }}{{mathrm{2}}}{mathrm{ + }}frac{{uppi }}{{mathrm{2}}}frac{{{T} – T_{mathrm{u}}}}{{T_{mathrm{c}} – T_{mathrm{u}}}}} right)} right)} hfill & {T_{mathrm{u}}{ < T} le T_{mathrm{c}}} hfill end{array}} right.,$$ (19) where FU8 is the rate of forcing for Model F8, T is the daily mean temperature (°C), TL = 4, Tu = 25, and Tc = 36. Analysis We assessed the ability of each chilling model to represent long-term trends in the chilling conditions by calculating CA using each chilling model for each station for 1951–2018. CA was calculated as the sum of CUi from 1 November in the previous year to 30 April. The trend of CA at each station was visualized as the slope of the linear regression of CA against year. We also calculated Pearson’s r between each pair of chilling models for each station to determine if the chilling models were interrelated. We calculated HR and CA of spring events for each species, station, and year to determine if the chilling models are consistent with the physiological assumption that the reduction in chilling would increase HR (Fig. 1). HR was calculated as the sum of FUi from 1 January to the date of onset of spring events using Model F1, and the performances of the other forcing models (F2–F8) were also tested. We also compared 1 January with the other two starting dates of temperature accumulation (15 January and 1 February) to test any potential difference causing by the date when temperature accumulation begins. CA was calculated as the sum of CUi from 1 November in the previous year to the date of onset of spring events. We chose 1 November as the start date for CA because the endodormancy of temperate trees began around 1 November52. We only tested the linear relationship because the data were better fitted by the linear regression than the exponential model (Fig. 3), even though CA was linearly or nonlinearly negatively correlated with HR17. Pearson’s r between CA and HR for all records was calculated for each species, with a significantly negative Pearson’s r (p  More

Acoustic data and analysis
Audio data were collected from online sound archives (Xeno-Canto—https://www.xeno-canto.org—and Macaulay libraries—https://www.macaulaylibrary.org), creating a pool of over 2000 audio tracks. We assessed the sound quality of these audio tracks by listening and through visual inspection of sound spectrograms. To capture intra-specific variation, we limited audio extraction to one drum per audio track (which also avoided pseudoreplication) and only included species for which at least 3 high-quality drums could be extracted. We retained 736 high-quality drums suitable for further analyses. These drums were distributed among 92 species (out of the 217 recognized species of woodpeckers50 and 22 genera, providing a representative sampling of the phylogenetic diversity found in this family (Fig.1b)). Background noise and other artifacts were reduced by wavelet continuous reconstruction (R ‘WaveletComp’ package72), following the methods and description outlined in previous work73. The full script is available on demand. Finally, 22 acoustic variables were extracted from these filtered sound samples using the R ‘Seewave’ package74. Given the pulsed-like nature of drumming, the chosen variables emphasized the temporal and amplitude-related (all normalized to the maximal amplitude within a given drumming signal) features of the sounds (Supplementary Table 1). These 22 variables were z-scored and then used in all subsequent analyses. Since these variables were partly correlated to varying degrees, we performed a principal component analysis (PCA) to reduce the number of descriptive variables quantifying drumming acoustic structure. This dimensionality reduction was useful for visualization and necessary for regularization (i.e. to prevent overfitting) in many of our analyses. This resulted in six principal components (PCs) with eigenvalues >1 which together explained 75% of the variance (Supplementary Table 11).
We used these variables to evaluate the similarity between species-specific drums, by performing a hierarchical cluster analysis (HCA)75,76 based on Euclidean distances and the ‘Ward.D2’ method (‘NbClust’ R package77). NbClust provides a clustering output resulting from the use of multiple indices (in the case of our analysis, 26 indices were used). The best number of clusters is chosen according to the majority rule, i.e. it is the one supported by the highest number of indices used. This entails creating a vector of acoustic features (22 raw acoustic measures or 22 PCs) for each of the 92 species in our dataset and calculating the Euclidean distance between these vectors to evaluate how close acoustically species were. Note that one can still use Euclidian distances in a non-orthonormal space to calculate the ‘distance’ between signals. The result is a distance metric that might give more weights to measures that co-vary. This could theoretically affect the clustering results. However, when we performed the same analysis using the 22 PCs, we obtained the same grouping (6 clusters) with very minor differences in species grouping and distances between clusters as shown in the relative length of branches (Supplementary Fig. 13). The output of this HCA established an optimal classification of woodpeckers’ drums into 6 main drumming types (Fig. 2a), described as follows:
– Acceleration (AC): Beak strikes decrease in amplitude as they are produced within successively shorter time intervals.
– Regular sequence (RS): Beak strikes are produced in bouts, each comprising a relatively fixed (stereotyped) number of strikes.
– Irregular sequence (IS): Beak strikes are produced in bouts, each comprising a variable number of strikes (as opposed to RS).
– Steady fast (SF): Beak strikes are produced with constant time intervals and at a similar amplitude, with a high pulse rate (on average >20 strikes/s).
– Steady slow (SS): Beak strikes are produced with constant time intervals and at a similar amplitude, with a low pulse rate (on average 2; see models Supplementary Tables 3–6), a likelihood ratio test (LRT) was conducted to test for the specific effect of predictor variables (Supplementary Table 4). Because both models (null and fitted) differ in their fixed effects, model comparison was performed on models fit by maximum likelihood (ML) with the phylogenetic correlation structure (Pagel’s λ) fixed to the estimates obtained from initial fit by REML. The statistics reported for model comparison are likelihood ratios.
PGLS models testing for a relationship between life-history variables and drumming structure included either of PC1 to PC6 as the dependent variable to investigate whether differences exist between these proxies for acoustic structure. Similarly, LDs were used to verify our results with these different loading combinations of drumming acoustic variables. No significant correlations were found between life-history traits and acoustic structure using LDs instead of PCs (Supplementary Table 6), indicating that the combination of structural variation captured by the LDs differed from that of the PCs, while not leading to fundamentally different conclusions. Similarly, no significant correlations were found between life-history traits and information content (no decrease in AICc >2; Supplementary Table 3), overall emphasizing that none of the variables investigated here (and which could have potentially affected drumming structure) seemed to have influenced species-specific information, or at least not directly.
Ancestral states reconstructions
We carried out two types of ancestral state reconstructions: discrete reconstruction of drumming status in Fig. 1b and of drumming types in Fig. 3a (using ‘ace’ from the R ‘ape’ package81), or continuous character reconstruction of drumming acoustic structure based on Brownian motion models (using ‘fastanc’ from the R ‘phytools’ package82 in Supplementary Figs. 4 and 5) and using relaxed Brownian motion model (using ‘rjmcmc.bm’ from the R ‘geiger’83 package in Fig. 4a and Supplementary Figs. 8 and 9).
While evaluating the likelihood that drumming was already present at an early stage of woodpecker’s phylogeny, we tried to represent the most complete tree of the family, based on very recent molecular data50. Note that strictly speaking, we evaluate the state at the root but at the next internal node, i.e. at the node including Picumninae and Picinae (the largest pie-chart in our Fig. 1b), as Wrynecks do no drum, and neither do honeyguides or barbets). To include species with an unknown drumming status in this discrete reconstruction, we attributed equal probability distribution between the 3 states (i.e. when the ‘drummer state’ of a species is unknown, the species is given, prior to ancestral state reconstruction, a 1/3 probability of belonging to each of the three categories ‘drummer’, ‘occasional drummer’ and ‘non-drummer’). Stochastic mapping was performed under an MCMC model, sampling the rate matrix from its posterior distribution for Q (‘Q = mcmc’ in make.simmap function from the R ‘phytools’ package), with an equiprobable default prior at the root, and 200 simulations. Under a symmetrical model for the probability to change among the three states, scaled likelihood on woodpeckers’ ancestral node indicated 56.4%, 38.3% and 5.3% probabilities of being a drummer, an occasional drummer and a non-drummer, respectively. This is in line with the fact that morphological adaptations for drilling (including reinforced rhamphotheca, frontal overhang and processus dorsalis pterygoidei) evolved in the ancestral lineage of Picumninae and Picinae64.
To prevent overfitting, the discrete reconstructions for drumming types were estimated for six different rate models: equal rate model (ER), symmetric rate model (SYM), all rates difference model (ARD) and three sequential transition models based on the normalized MIL as measures of complexity as shown in Supplementary Fig. 7 ((SF leftrightarrow SS leftrightarrow DK leftrightarrow AC leftrightarrow RS leftrightarrow IS)). These three models assumed (1) sequential and equal, (2) sequential and incremental and (3) sequential and reversed transition rates, respectively. The number of parameters for these 6 rate models were 25, 1, 15, 1, 2 and 10. The final regularized likelihoods of each ancestral states were then obtained by model averaging using Akaike weights.
Calculation of information at different evolutionary steps was carried out as an extension of the drumming types reconstruction described above. From the discrete ancestral reconstruction procedure, probability distributions of drumming types were obtained for each node of the phylogenetic tree. We then obtained probability distributions at 20 fixed time intervals (dt = 1 myr) by linear interpolation. Using these probability distributions, we sampled drumming types proportionally from extant species descending the node closest to the time interval to estimate ancestral information values. This bootstrap procedure was repeated 30 times in order to obtain reliable estimates of mean and standard error. In this manner, we obtained information-through-time plots. These plots quantify a putative diversity of drumming signals in the clade at a particular point in time. They are similar in spirit to the disparity-through-time plots that have been used to measure specific morphological diversity in a clade through time using phylogenetic trees based on molecular data in combination with morphological measures in extant species84.
Continuous ancestral character trait reconstruction of drumming acoustic structure was carried out using either the six PCs that explain variation among the 22 drumming acoustic variables, or the six LDs that explain the variation in discriminating potential among the same variables (see above, ‘Acoustic data and analysis’ and ‘Calculation of information’ sections; Supplementary Figs. 4 and 5). The results and conclusions were similar for all PC’s and since the PC1 component has strong loading of multiple acoustic variables and the highest acoustic structure variance explained (Supplementary Table 11) it serves well as an illustrative example. The measure of phylogenetic signal on continuous traits (i.e. the historical contingency between species-specific drums that renders a trait non-randomly distributed along the phylogenetic tree) was made using Pagel’s lamba (Supplementary Table 2).
Reconstructing information content from raw MIL values would not have been biologically relevant since information calculation is based on the number of species involved, a factor that changes as branches merge going backward along the phylogenetic tree. We thus reconstructed MIL based on the normalized MIL values to avoid this pitfall. We used a Bayesian model implemented in the R package ‘Geiger’83 (model ‘rbm’ in the function ‘rjmcmc.bm’) to estimate branch-specific rates of trait evolution (i.e. changes in rates through time and across lineages). In this method, a reversible jump Markov Chain Monte Carlo (MCMC) sampling algorithm is used to detect shifts in rates of continuous traits evolution under a relaxed Brownian motion model85. The results of the model fit were summarized by the branch-specific average rate, estimated from the posterior samples. To obtain relative variations in posterior average rates, drumming structure (PC1-PC6) and MIL were standardized, i.e. these traits were divided by their standard deviation prior to running the ‘rbm’ models.
Analytical simulations of selection for information
In Fig. 3c, we compared that reconstructed evolution of information to what might be expected in different scenarios to further support those conclusions. More specifically, we estimated the ancestral MI for two simulated scenarios using an analytical model that describes species-specific information based on the probability of correct detection and the number of species (see ‘Calculation of information’ section). In the ‘No Diversifying Selection’ scenario (dark brown), the probability of correct detection for the initial pair of species, p2, is first estimated from the data using the approach described in the main text. It is then assumed that that additional species are randomly just as different/similar than these original species pair, yielding a probability of correct detection through time given by (p_c(t) = p_2^{n_s(t) – 1}), where ns(t) is the number of species at a given time. In the ‘Strong Diversifying Selection’ scenario (light brown), the probability of correct detection estimated at the first time point in our reconstruction (−20 M years ago, 3 species) is kept constant, (p_c(t) = p_2). In other words, the only species that survive would be species that can discriminate themselves from all other species equally well than the currently existing species. The reconstructed (actual) scenario is found between these two extreme values, showing that the drumming types are clearly not random but were also not under high evolutionary pressure to increase species-specific information. New drumming types evolved and species within types used signals that were distinct enough to result in the maintenance of normalized MI.
In Supplementary Fig. 6, we showed that the non-normalized reconstructed MI increased more rapidly when new drumming types appeared but that the normalized MI was relatively constant, reflecting the fact that the appearance of novel drumming types could co-occur with rapid radiation and increase in species numbers.
Playback experiments
Initial preparation involved identifying and mapping the areas prone to high densities of great-spotted woodpeckers Dendrocopos major, the study species of this experimental phase, using GIS maps provided by the LPO (French Bird Protection Organization). D. major is commonly found in European forests, ranging from open coniferous to mature deciduous forests. Playback experiments were carried out on wild individuals around Saint-Etienne, France, during this species’ breeding season (February–April 2017). All experiments were performed in accordance with relevant guidelines and regulations including French national guidelines, permits and regulations regarding animal care and experimental use (approval no. D42-218-0901, ENES lab agreement, Direction Départementale de la Protection des Populations, Préfecture du Rhône).
Two sets of experiments were conducted over the course of the breeding season, although we implemented the same general design which consisted in simulating a territorial intrusion. Playback stimuli tracks consisted of eight drums spread unevenly over about 60 s, aiming at representing the variation encountered in natural sequences (ref. 44 and personal observations). The first experiment (Exp. 1) aimed at investigating D. major’s response to conspecific vs. heterospecific drums. The other experiment (Exp. 2) aimed at investigating D. major’s response to drums from conspecifics vs. drums modified through acoustic manipulation (i.e. signal re-synthesis). D. major typically drums with an ‘acceleration’ pattern, which is mainly characterized by a shortening of the inter-strike time interval, a progressive decrease in strikes’ amplitude, and a gradual change in spectral properties as strikes get faster and weaker.
In Exp. 1, we used a paired and randomized order design, presenting each focal individual with one D. major drum and one drum from one out of 4 different species: 2 of which have very different drumming patterns (Picus canus and Dryobates minor, both producing ‘steady fast’ drums), and 2 others which have similar (accelerating) drumming patterns (Dendrocopos syriacus and Dendrocopos hyperythrus). A potentially confounding factor (which is nevertheless in line with our phylogenetic analyses) lies in that the allopatric species producing a drum similar to that of D. major also happened to be closely related to our model species. We carried out 48 playback experiments (testing 24 individuals with one of 4 categories of paired signals).
In Exp. 2, we altered one of the 3 acoustic features described above or all of them together (thus having 4 categories of modified signals), using Praat sound analysis software86. The design was paired so that each focal individual was exposed to one conspecific drum and one modified drum, following a randomized presentation order. This led to 48 playback experiments (24 individuals, each tested with one of 4 categories of paired signals).
Within each of Exp. 1 and Exp. 2, tested individuals were all separated by at least 500 m, ensuring different identities since their territory sizes vary between 200 and 400 m87,88. Upon visual or aural detection of (an) active individual(s), the experimenter set up an Anchor Megavox loudspeaker at about 1–1.5 m from ground level. The speaker was connected to an Edirol R-09 recorder (stimuli tracks were created and stored as WAV files, 44.1 kHz sampling frequency). Playbacks started at about 50 m from where the experimenter last saw or heard the focal individual. Following the work from Schuppe and colleagues89, playback intensity was calibrated and kept at about 80 dB measured 1 m away from the speaker. Behavioural data collection started when the first drum of the stimuli track was broadcasted and lasted 10 min from that moment on. To document focal individuals’ responses, notes were taken manually and continuously, while audio was recorded with a Sennheiser ME67 microphone mounted on a tripod and connected to a digital recorder (Zoom H4N, 44.1 kHz, 16 bit). If a response was elicited from multiple individuals in the area, only the one from a particular individual (ideally the one seen or heard before setting up the experiment) was monitored and used in further analyses. Six behavioural variables were reported, namely the number of screams, the number of drums, the approach (which was divided into three categories: ‘within 25 m’, ‘25–50 m’ and ‘further than 50 m’) as well as the latencies to first scream and drum and the latency to closest approach. When no occurrence was observed for the first three behaviours, latencies were set by default to the maximum value, i.e. the duration of the full experiment (10 min = 600 s). To characterize D. major’s behaviour, a PCA was then performed on scaled/centred data, where we retained the first principal component (‘Playback-PC1’) as an indicator of the behavioural response’s strength. A higher Playback-PC1 score indicates a stronger territorial response, i.e. more screams, a closer approach to the speaker and shorter latencies to these 2 behaviours. A second significant component resulted from this PCA, which represented the drumming’s response (inversely related: a higher Playback-PC2 score indicates fewer drums and a longer latency to drum; see Supplementary Table 13). None of the pairwise comparisons were statistically significant for PC2, besides a stronger drumming response to drums resynthesized without temporal variation than to D. major drums (Supplementary Fig. 15). This can be explained by the fact that birds were tested during their breeding season. At this time, drumming behaviour is likely to occur more consistently and commonly across experiments, independently from the stimulus played back, while screams and approach do not occur unless threat of intrusion is clear. Therefore, we used Playback-PC1 to represent birds’ behavioural response in our analysis (as it is in addition explaining much more variance in the behavioural data than Playback-PC2). Note that, as two playback sets were involved in this study, while we considered them independently in our statistical analysis, for standardization of the behavioural scale, we used the same polynomial equation. More specifically, the linear equation obtained from the loading scores of Exp. 1 was applied to the behavioural data of Exp. 2 for computation of Playback-PC2 scores.
Finally, distances were approximated during continuous note-taking and confirmed post-experimentally using a National Geographic 4*21 rangefinder (measurement accuracy: ±1 m up to 200 m). Sex was not documented as sometimes birds were not seen (but just heard drumming or calling back at our playback), which we nevertheless believe to be negligible since both sexes drum and are territorial in this monogamous species44,90.
Statistical analyses tested for differential responses of focal birds to drums of their own species versus either another species or a modified resynthesized condition. A paired comparison design was used by means of LMMs and contrasts using R software (‘lme4’ and ‘lsmeans’ packages)91,92. LMMs included study day and time, order of presentation and focal bird identity as random factors, and tested for a fixed effect of the interaction between treatment and group of paired condition. Contrasts were then computed between treatments (i.e. conspecific versus non-conspecific drums) for each group (i.e. each paired testing condition, such as D. major versus D. minor for which n = 6 birds were exposed to paired playback presentations—see Fig. 5a, b). Before contrasts and using the ‘lsmeans’ function, a Tukey adjustment for multiple testing was used; two-sided statistics are reported.
Reporting summary
Further information on research design is available in the Nature Research Reporting Summary linked to this article. More

Publication collection and data compilation
The detailed methods of publication search and data collection were described in our related work7. In brief, 10 databases in Web of Science (WoS; 1 January 1900 to 13 December 2016) including BIOSIS Previews, Chinese Science Citation Database, Data Citation Index, Derwent Innovations Index, Inspec, KCI-Korean Journal Database, MEDLINE, Russian Science Citation Index, SciELO Citation Index, and WoS Core Collection were used for searching peer-reviewed publications that reported GCMEs. The 18 keywords for WoS title search were: global change, climate change, free-air carbon dioxide enrichment, free-air CO2 enrichment, elevated carbon dioxide, elevated CO2, elevated atmospheric CO2, CO2 enrichment, eCO2, [CO2], warming, elevated temperature, changing precipitation, increased precipitation, decreased precipitation, nitrogen deposition, nitrogen addition, and nitrogen application. Through these search, 310,177 publication records that might be relevant to our topic were found.
First, we identified all the 310,177 records via reading each title. Second, we read the abstracts of all the records collected in the first step to further screen publications. During the two steps, we excluded 291,436 records because these studies were reviews/meta-analyses or conducted in non-terrestrial ecosystems such as oceans. Third, we read the methods of the remaining 18,741 publications to identify which of them met the following three inclusion criteria:
1.
Publications reported results of outdoor GCMEs which had at least three control and global change treatment plots ( > = 1 m2).

2.
The GCMEs were conducted in terrestrial ecosystems except for croplands and lab incubation studies.

3.
The GCMEs aimed to examine effects of simulated global change drivers on carbon, nitrogen, and water-cycle variables as well as plant and microbial parameters.

During the screening in the third step, 1,290 publications met these defined criteria.
We subsequently cross-checked the list of the 1,290 publications with references cited by the previous review/meta-analysis articles in global change research as well as the 1,290 publications, and collected 756 publications. In addition, 184 studies were collected by searching the websites of ecology laboratories and experiment networks and checking the references of the papers downloaded from these websites. In total, 2,230 publications were collected in the original version of the database7. Moreover, another 12 publications were found when we checked and reorganized all the data extracted from the 2,230 publications8,9,10,11,12,13,14,15,16,17,18,19. This database compiled 11 plant production and ecosystem carbon exchange variables including net primary productivity (NPP), above- and below-ground NPP (ANPP and BNPP), total biomass, aboveground biomass (AGB), root biomass, litter mass, gross and net ecosystem productivity (GEP and NEP), and ecosystem and soil respiration (ER and SR). Data of mean values, standard deviations or standard errors, and sample sizes (number of plot replications) of these variables in the control and treatment (e.g., elevated CO2, nitrogen addition, warming, increased/decreased precipitation, or their combinations) groups were extracted from each publication when possible. The figures were digitized using SigmaScan Pro 5.0 (SPSS, Inc.) and the numerical values were extracted when a publication presented experimental data graphically. Data of the experiments that were conducted over less than one year/growing season were excluded in this database. However, we included short-term data from tundra studies because most of measurements in this ecosystem were performed during July-August. Overall, 5,213 pairs (the control versus global change treatment) of plant production and ecosystem carbon exchange samples were collected in this database, having 2,247, 2,120, 81, and 765 pairs from single-, two-, three-, and four-factor manipulative experiments, respectively (Fig. 1).
Fig. 1

Number of samplings. Number of sample pairs of ecosystem carbon-cycling variables including net primary productivity (NPP), above- and below-ground NPP (ANPP and BNPP), total biomass, aboveground biomass (AGB), root biomass, litter mass, gross and net ecosystem productivity (GEP and NEP), and ecosystem and soil respiration (ER and SR) extracted from publications reporting single-, two-, three-, and four-factor global change manipulative experiments.

Full size image

Environmental metadata: Climate and vegetation
Information on the locations and altitudes of each experimental site, site climate including mean annual temperature (MAT) and precipitation (MAP) as well as wetness index ((frac{{rm{MAP}}}{{rm{MAT+10}}})), ref. 20, and vegetation types were extracted from each of the 2,242 publications. If a study did not report climate characteristics for its experimental site, data of MAT and MAP were downloaded from Climate Model Intercomparison Project phase 5 (CMIP5; https://esgf-node.llnl.gov/projects/cmip5/) based on the site coordinate. The dataset selection in CMIP5 was “historical (simulation of recent past 1850–2005)” and the climate data averaged from 20 (i.e. BCC_CSM1_1, BCC_CSM1_1_M, CANESM2, CCSM4, CMCC_CM, CMCC-CMS, CNRM-CM5, CSIRO_MK3_6_0, GFDL_CM3, GISS_E2_H, HADGEM2_AO, HADGEM2_ES, INMCM4, MIROC_ESM, MIROC_ESM_CHEM, MIROC5, MPI_ESM_LR, MPI_ESM_MR, MRI_ESM1, and NORESM1_M)21, that contained historical climate data, out of the 35 global climate models available in CMIP5 were used in this study. In addition, we downloaded data of climate means at global 1 × 1° land grid cells from Princeton University (http://hydrology.princeton.edu/data/pgf/v3/) to construct global climate space. Moreover, we classified ecosystems subjected to ecosystem manipulative experiments into five typical types: forests (mature forests and tree seedlings), grasslands (grasslands, meadows, short- and tall-grass prairies, temperate/semi-arid steppes, shrublands, savannas, pastures, and old-fields), tundra, wetlands (peatlands, bogs, marshes, and fens), and deserts.
Metadata of experimental facilities and performance
Information on CO2 enrichment and warming facilities were also extracted from the related publications reporting CO2 and warming effects on plant production and ecosystem carbon exchange. Facilities used in elevated CO2 experiments included greenhouse, open-top chamber, free-air CO2 enrichment, and tunnels. Warming experiments primarily used greenhouse, open-top chamber, soil heating cables, infrared radiator, and infrared reflector to elevate vegetation canopy and soil temperature. In addition, the manipulation magnitudes of global change drivers imposed by manipulative experiments, such as the increases in CO2 concentrations (ppm) and temperature (°C), the changes in precipitation amount (mm), and the rates of nitrogen input (g N m−2 yr−1), were also collected and added into this updated database. More

The decision by an individual to move from one area to another is mediated by a number of factors, such as resource quality and availability, predation risk and local environmental conditions, all of which will influence its survival and reproductive output1,4. The challenge for conservationists is understanding how these individual decisions can affect population dynamics, home ranges and ultimately species’ survival1.
Home ranges of carnivores should overlap and in some cases envelop those of their prey species. Southern Ground-hornbills feed on a variety of prey, ranging from snakes, rabbits and birds to invertebrates12,18. Through tracking Southern Ground-hornbill movements, we were able to show that group home ranges during the early and late dry seasons were larger than in the wet season. As the Southern Ground-hornbill breeding season in South Africa coincides with the warm, wet summer months, prey availability, especially that of invertebrates, is expected to be higher20,21, suggesting that individuals would not need to travel as extensively to find sufficient food. Furthermore, in the late dry season, groups used between 76 and 115% of their home ranges. This was likely a result of having to increase their search for food and relaxation of the central place foraging required around the nest during the breeding season.
Previous research on Southern Ground-hornbill home ranges has recorded group densities ranging from one group per 4000 ha (communal areas in Zimbabwe22), to one group every 10,000 ha (KNP14), with one group in the Limpopo Valley having a home range close to 20,000 ha21. These results were obtained by direct observations of active nest sites or using VHF radio transmitters. In our study using GPS data, we showed that home range sizes of Southern Ground-hornbills within KNP vary considerably. Despite this, our results confirmed the findings of Theron et al.21 and Zoghby et al.20, demonstrating a restricted and contracted home range during the breeding season, when group movements are concentrated around the nest site (central place foraging). Presumably, breeding success would influence the extent of wet seasonal home range for Southern Ground-hornbills, with groups abandoning their central place foraging behaviour when nests fail. Wyness19 reported that of four Southern Ground-hornbill groups studied in the Association of Private Nature Reserves (APNR) adjacent to the KNP, the three that bred successfully in the year of their study showed a breeding season range reduction to between 24–36% of their non-breeding home range. The unsuccessful group used 70% of their home range during this time19. Surprisingly, the groups within the KNP did not show such a definitive pattern in home range size reduction associated with breeding success, although all groups that attempted breeding did show a wet seasonal home range reduction. Of the six Southern Ground-hornbill groups monitored in our study, four groups bred successfully, one group’s attempt failed (Ngotso Camp), and the breeding status for the third group (Shingwedzi) was unknown. The groups that bred successfully used 21–97% of their respective home ranges, with the unsuccessful group using 85% of their home range (See Table 1).
Southern Ground-hornbills are known to favour more open habitats for foraging20,23. Our results supported this, with groups selecting the open woodland and grassland habitat types year-round, following their availability within the landscape.
Although Southern Ground-hornbill seasonal territory size differed significantly amongst the groups, they all showed a decrease in the amount of low shrubland and an increase in the amount of grassland habitat used with increased territory size. Similarly, as seasonal territory sizes increased, the amount of low-medium woody cover (25–50%) decreased. Thus, when selecting an area for a reintroduction of Southern Ground-hornbill groups, the ratio of low-medium woody cover (low shrubland) to grassland, calculated based on the national land cover datasets available, should be taken into account, as this will likely influence the home range size and the number of groups that could be supported in an area.
Although an understanding of the changes and restrictions in territory size is important for the management of a species, the types of movements adopted within a population will influence the management actions needed for their conservation, such as ensuring connectivity or access to certain resources1. Conservation policy and management actions are less effective when interventions do not integrate both the spatial and temporal changes in habitat use and the scale of species movements1,3. The results from the first-passage time analysis of Southern Ground-hornbill movements showed that the different groups did not consistently demonstrate seasonal patterns in the scale at which they concentrated their foraging efforts. The mean distances travelled for all trajectory paths, classified as active foraging behaviour, were similar and lower in the late wet and early dry seasons compared with the late dry and early wet seasons. Movement between foraging resource patches or mean relocation distances were highest in the wet season months, with the maximum mean distances travelled during the early wet season and the start of the breeding period. Overall prey abundance for Southern Ground-hornbills is generally higher in the wetter months, resulting in a decrease in relocation distances. Our results support the theory that Southern Ground-hornbill wet season movements are most likely influenced by the need to travel to and from the nest site to provision prey to the incubating female and growing nestling. Once resources closer to the nest are depleted, the distances travelled to access additional habitats and prey would likely increase.
Southern Ground-hornbills seemingly prefer nest sites surrounded by more open woodland habitat24,25. Habitat structure and the diversity of habitat types within a 3 km radius around the nest site positively influenced Southern Ground-hornbill nesting success. An increase in the density of woody habitat surrounding the nest site, however, had a negative impact on Southern Ground-hornbill breeding success24, possibly owing to decreased foraging opportunities, an increased risk of predation or an increase in foraging effort beyond a value which is beneficial.
Habitat structure will likely promote or inhibit the types of movement that can occur in an area. The results from the multinomial regression (Table 5) indicate that the likelihood of a movement behaviour being classified as “foraging” within the open woodland, grassland and dense thicket habitat types was higher than the behaviour being attributed to “relocating”. This is to be expected for open woodland, and grassland habitats as these are both ideal open foraging habitats for Southern Ground-hornbills20,23 and are used year-round in proportion to their availability. Southern Ground-hornbills spend around 70% of their day walking12 and have been shown to travel distances of up to 10.6 km in a day20. Having to navigate through dense thicket vegetation in an area may increase the amount of time spent there, possibly accounting for why this habitat type is predicted to be used more for “foraging”-type behaviour as opposed to “relocating” behaviour. Travel through areas of low shrubland habitat was considered “relocating” behaviour, suggesting that within this habitat type, it is more profitable for Southern Ground-hornbills to move further, and the corresponding chance of finding food greater, than conducting area-restricted searches and spending longer periods concentrated in one patch.
When comparing movements between habitats allocated to “resting” as opposed to “foraging”, the time spent in all habitats was most likely as a result of “foraging”. As GPS locations were only recorded during the day, switching off at dusk (~ 18h00) when Southern Ground-hornbills would roost for the night, habitat preferences for “resting” movements may not have been recorded. Moreover, during the day, Southern Ground-hornbills may not be actively selecting for specific habitat types in which to roost or rest. They may simply be roosting or resting at a chosen site to escape the midday heat within the habitat type in which they were “foraging” or “relocating”.
We were unable to explore differences in movement relating to specific characteristics of the tagged bird (age, sex, helper versus breeder status, etc.) in our study. However, future research should consider study designs able to account for these potential differences, as García-Jiménez et al.26 showed that both the breeding season and sex of the individual influence displacement and distance travelled in Pyrenean Bearded Vultures (Gypaetus barbatus). They found that all individuals travelled more in the breeding season, with females having greater cumulative and maximum distances regardless of the season. More

1.
Lomolino, M. V., Riddle, B. R. & Brown, J. H. Biogeography (Sinauer Associates Inc, Sunderland, 2006).
Google Scholar 
2.
Chamorro, D., Olivero, J., Real, R. & Muñoz, A.-R. Environmental factors determining the establishment of the African Long-legged Buzzard Buteo rufinus cirtensis in Western Europe. Ibis (Lond. 1859). 159, 331–342 (2017).
Google Scholar 

3.
Pearson, R. G. & Dawson, T. P. Predicting the impacts of climate change on the distribution of species: Are bioclimate envelope models useful ?. Glob. Ecol. Biogeogr. 12, 361–371 (2003).
Google Scholar 

4.
Sun, J. et al. Modeling the potential distribution of Zelkova schneideriana under different human activity intensities and climate change patterns in China. Glob. Ecol. Conserv. 21, e00840 (2020).
Google Scholar 

5.
Mackey, B. & Lindemayer, D. Towards a hierarchical the spatial distribution of animals. J. Biogeogr. 28, 1147–1166 (2011).
Google Scholar 

6.
Márquez, A. L., Real, R., Olivero, J. & Estrada, A. Combining climate with other influential factors for modelling the impact of climate change on species distribution. Clim. Change 108, 135–157 (2011).
ADS  Google Scholar 

7.
IPCC. Climate Change 2014: Synthesis Report. Contribution of Working Groups I, II and III to the Fifth Assessment Report of the Intergovernmental Panel on Climate Change (2014) https://doi.org/10.1017/CBO9781107415324.004.

8.
Root, T. L. et al. Fingerprints of global warming on wild animals and plants. Nature 421, 57–60 (2003).
ADS  CAS  PubMed  Google Scholar 

9.
Karl, T. R. & Trenberth, K. E. Modern global climate change. Science (80-). 302, 1719–1723 (2003).
ADS  CAS  Google Scholar 

10.
Walther, G.-R. et al. Ecological responses to recent climate change. Nature 416, 389–395 (2002).
ADS  CAS  PubMed  Google Scholar 

11.
Lazo-Cancino, D. et al. The impacts of climate change on the habitat distribution of the vulnerable Patagonian-Fueguian species Ctenomys magellanicus (Rodentia, Ctenomyidae). J. Arid Environ. 173, 104016 (2019).
Google Scholar 

12.
Borzée, A. et al. Climate change-based models predict range shifts in the distribution of the only Asian plethodontid salamander: Karsenia koreana. Sci. Rep. 9, 11838 (2019).
ADS  PubMed  PubMed Central  Google Scholar 

13.
Thomas, C. D. & Lennon, J. J. Birds extend their ranges northwards. Nature 399, 213 (1999).
ADS  CAS  Google Scholar 

14.
Wallingford, P. D. et al. Adjusting the lens of invasion biology to focus on the impacts of climate-driven range shifts. Nat. Clim. Change. 10, 398–405 (2020).
ADS  Google Scholar 

15.
Castro, A., Muñoz, A.-R. & Real, R. Modelling the spatial distribution of the Tengmalm’s owl Aegolius funereus in its Southwestern Palaeartic limit (NE Spain). Ardeola 55, 71–85 (2008).
Google Scholar 

16.
Thomas, C. D. et al. Ecological and evolutionary processes at expanding range margins. Nature 411, 577–581 (2001).
ADS  CAS  PubMed  Google Scholar 

17.
Huntley, B. et al. Potential impacts of climatic change upon geographical distribution of birds. Ibis (Lond. 1859). 148, 8–28 (2006).
Google Scholar 

18.
Maclean, I. M. D. et al. Climate change causes rapid changes in the distribution and site abundance of birds in winter. Glob. Change. Biol. 14, 2489–2500 (2008).
Google Scholar 

19.
Elorriaga, J. & Muñoz, A.-R. First breeding record of North African Long-legged Buzzard Buteo rufinus cirtensis in continental Europe. Br. Birds 103, 399–401 (2010).
Google Scholar 

20.
Nieto, I., Chamorro, D., Palomo, L. J., Real, R. & Muñoz, A.-R. Is the Eurasian Reed Warbler a regular wintering species in the Iberian Peninsula? Ringing data say yes. Acta Ornithol. 53, 61–68 (2018).
Google Scholar 

21.
Ambrosini, R. et al. Climate change and the long-term northward shift in the African wintering range of the barn swallow Hirundo rustica. Clim. Res. 49, 131–141 (2011).
Google Scholar 

22.
Møller, A. P., Rubolini, D. & Lehikoinen, E. Populations of migratory bird species that did not show a phenological response to climate change are declining. Proc. Natl. Acad. Sci. U. S. A. 105, 16195–16200 (2008).
ADS  PubMed  PubMed Central  Google Scholar 

23.
Real, R., Márquez, A. L., Olivero, J. & Estrada, A. Species distribution models in climate change scenarios are still not useful for informing policy planning: An uncertainty assessment using fuzzy logic. Ecography (Cop.) 33, 304–314 (2010).
Google Scholar 

24.
Root, T. L. & Schneider, S. H. Conservation and climate change: The challenges ahead. Conserv. Biol. 20, 706–708 (2006).
PubMed  Google Scholar 

25.
Hahn, S., Bauer, S. & Liechti, F. The natural link between Europe and Africa—2.1 billion birds on migration. Oikos 118, 624–626 (2009).
Google Scholar 

26.
Evans, P. R. & Lathbury, G. W. Raptor migration across the Strait of Gibraltar. Ibis (Lond. 1859). 115, 572–585 (1972).
Google Scholar 

27.
Bijlsma, R. G. Bottleneck Areas for Migratory Birds in the Mediterranean Region: An Assessment of the Problems and Recommendations for Action. (ICBP Study report no. 18, 1987).

28.
Gantenbein, B. & Largiadèr, C. R. The phylogeographic importance of the Strait of Gibraltar as a gene flow barrier in terrestrial arthropods: A case study with the scorpion Buthus occitanus as model organism. Mol. Phylogenet. Evol. 28, 119–130 (2003).
CAS  PubMed  Google Scholar 

29.
García-Mudarra, J. L., Ibáñez, C. & Juste, J. The Straits of Gibraltar: Barrier or bridge to Ibero-Moroccan bat diversity?. Biol. J. Linn. Soc. 96, 434–450 (2009).
Google Scholar 

30.
Carranza, S., Harris, D. J., Arnold, E. N., Batista, V. & de la Gonzalez, V. Phylogeography of the lacertid lizard, Psammodromus algirus, in Iberia and across the Strait of Gibraltar. J. Biogeogr. 33, 1279–1288 (2006).
Google Scholar 

31.
de Juana, E. & Comite de Rarezas de la Sociedad Española de Ornitología, (SEO). Observaciones de Aves raras en España, año 1995. Ardeola44, 119–141 (1997).

32.
de Juana, E. & Comité Ibérico de Rarezas de la Sociedad Española de Ornitología, (SEO). Observaciones homologadas de aves raras en España y Portugal. Informe de 1992. Ardeola41, 103–117 (1994).

33.
Copete, J. L. et al. Observaciones de Aves Raras en España, 2012 y 2013. Ardeola 62, 453–508 (2015).
Google Scholar 

34.
Dies, J. I. et al. Observaciones de Aves Raras en España, 2008. Ardeola 57, 481–516 (2010).
Google Scholar 

35.
Elorriaga, J. & Muñoz, A.-R. Hybridisation between the Common Buzzard Buteo buteo buteo and the North African race of Long-legged Buzzard Buteo rufinus cirtensis in the Strait of Gibraltar: Prelude or preclude to colonisation?. Ostrich J. Afr. Ornithol. 84, 41–45 (2013).
Google Scholar 

36.
Chamorro, D., Olivero, J., Real, R. & Muñoz, A.-R. El cambio del clima y la Barrera Biogeográfica del Estrecho de Gibraltar para las aves africanas. In Avances en biogeografía: Áreas de distribución: entre puentes y barreras (eds. Zotano, J. G., et al.) 172–180 (Tundra Ediciones, 2016).

37.
Ramírez, J. et al. Spring movements of Rüppell’s Vulture Gyps rueppellii across the Strait of Gibraltar. Ostrich 82, 71–73 (2011).
Google Scholar 

38.
del Hoyo, J., Elliott, A. & Sargatal, J. Handbook of the Birds of the World. Vol. 2. New World Vultures to Guineafowl. (Lynx Edicions, Barcelona, 1994).

39.
Ferguson-Lees, J. & Christie, D. A. Raptors of the World (Houghton Mifflin Harcourt, Boston, 2001).
Google Scholar 

40.
Irby, L. H. The Ornithology of the Straits of Gibraltar (Taylor & Francis, Abingdon, 1895).
Google Scholar 

41.
Ávila, A., Guirado, M. A. & Navarrete, J. Observaciones de Aves Raras. Busardo Moro. Ardeola 51, 548 (2004).
Google Scholar 

42.
Cramp, S. & Simmons, K. E. L. Handbook of the Birds of Europe, the Middle East and North Africa: The Birds of the Western Palearctic. Vol. 2: Hawks to Bustards. (Oxford University Press., Oxford, 1980).

43.
Zadeh, L. A. Fuzzy sets. Inf. Control 8, 338–353 (1965).
MATH  Google Scholar 

44.
Real, R. et al. Conservation biogeography of ecologically interacting species: The case of the Iberian lynx and the European rabbit. Divers. Distrib. 15, 390–400 (2009).
Google Scholar 

45.
Olivero, J., Real, R. & Márquez, A. L. Fuzzy chorotypes as a conceptual tool to improve insight into biogeographic patterns. Syst. Biol. 60, 645–660 (2011).
PubMed  Google Scholar 

46.
Olivero, J. et al. Recent loss of closed forests is associated with Ebola virus disease outbreaks. Sci. Rep. 7, 14291 (2017).
ADS  MathSciNet  PubMed  PubMed Central  Google Scholar 

47.
Robertson, M. P., Villet, M. H. & Palmer, A. R. A fuzzy classification technique for predicting species’ distributions: Applications using invasive alien plants and indigenous insects. Divers. Distrib. 10, 461–474 (2004).
Google Scholar 

48.
Salski, A. Ecological applications of fuzzy logic. In Ecological Informatics (ed. Recknagel, F.) 3–14 (Springer, New York, 2006).

49.
Real, R., Barbosa, A. M. & Bull, J. W. Species distributions, quantum theory, and the enhancement of biodiversity measures. Syst. Biol. 66, 453–462 (2017).
PubMed  Google Scholar 

50.
Acevedo, P. & Real, R. Favourability: Concept, distinctive characteristics and potential usefulness. Naturwissenschaften 99, 515–522 (2012).
ADS  CAS  PubMed  Google Scholar 

51.
Hosmer, D. W. & Lemeshow, S. Goodness of fit tests for the multiple logistic regression model. Commun. Stat. Theory Methods 9, 1043–1069 (1980).
MATH  Google Scholar 

52.
Thévenot, M., Vernon, R. & Bergier, P. The Birds of Morocco: An Annotated Checklist. BOU Checklist No. 20. (British Ornithologist’s Union, 2003).

53.
Muñoz, A.-R., Márquez, A. L. & Real, R. An approach to consider behavioral plasticity as a source of uncertainty when forecasting species’ response to climate change. Ecol. Evol. 5, 2359–2373 (2015).
PubMed  PubMed Central  Google Scholar 

54.
Muñoz, A.-R., Real, R., Barbosa, A. M. & Vargas, J. M. Modelling the distribution of Bonelli’s eagle in Spain: Implications for conservation planning. Divers. Distrib. 11, 477–486 (2005).
Google Scholar 

55.
Guisan, A. & Zimmermann, N. E. Predictive habitat distribution models in ecology. Ecol. Modell. 135, 147–186 (2000).
Google Scholar 

56.
Guisan, A. & Thuiller, W. Predicting species distribution: Offering more than simple habitat models. Ecol. Lett. 8, 993–1009 (2005).
Google Scholar 

57.
Gouveia, S. F. et al. Ecophysics reload—exploring applications of theoretical physics in macroecology. Ecol. Modell. 424, 109032 (2020).
Google Scholar 

58.
Cramp, S. Handbook of the Birds of Europe, the Middle East and North Africa: The birds of the Western Paleartic, Vol. VI. (Oxford University Press, Oxford, 1992).

59.
Hardin, G. The competitive exclusion principle. Science (80-). 131, 1292–1297 (1960).
ADS  CAS  Google Scholar 

60.
Jowers, M. J. et al. Unravelling population processes over the Late Pleistocene driving contemporary genetic divergence in Palearctic buzzards. Mol. Phylogenet. Evol. 134, 269–281 (2019).
PubMed  Google Scholar 

61.
Waters, J. M., Fraser, C. I. & Hewitt, G. M. Founder takes all: Density-dependent processes structure biodiversity. Trends Ecol. Evol. 28, 78–85 (2013).
PubMed  Google Scholar 

62.
Gil-Velasco, M. et al. Observaciones de Aves Raras en España, 2017. Ardeola 66, 169–204 (2019).
Google Scholar 

63.
Gutiérrez, R. et al. Observaciones De Aves Raras En España, 2011. Ardeola 60, 437–506 (2013).
Google Scholar 

64.
Gutiérrez, R. et al. Observaciones aves raras en España, 2010. Ardeola 59, 353–411 (2012).
Google Scholar 

65.
Jara, J. et al. Comité Português de Raridades. Anuário Ornitológico da SPEA 7, 140 (2010).
Google Scholar 

66.
Corso, A. Successful mixed breeding of Atlas Long-legged Buzzard and Common Buzzard on Pantelleria, Italy, in 2008. Dutch Bird. 31, 224–226 (2009).
Google Scholar 

67.
Väli, Ü et al. Widespread hybridization between the Greater Spotted Eagle Aquila clanga and the Lesser Spotted Eagle Aquila pomarina (Aves: Accipitriformes) in Europe. Biol. J. Linn. Soc. 100, 725–736 (2010).
Google Scholar 

68.
del Junco, O. & González, B. L. nueva especie de Vencejo en el Paleártico: Apus caffer. Ardeola 13, 115–127 (1969).
Google Scholar 

69.
Ferrero, J. J. Situación del Elanio azul (Elanus caeruleus) en el Mediterráneo. In Biología y Conservación de las Rapaces Mediterráneas,1994 (eds. Muntaner, J. & Mayol, J.) 101–115 (SEO, Sociedad Española de Ornitología, 1996).

70.
Ramírez, J., Simón, M., Solís, S., Pérez, C. & García, E. Vencejo Moro Apus affinis Observaciones de Aves Raras en España. Ardeola 49, 161 (2002).
Google Scholar 

71.
Molina, B., Prieta, J., Lorenzo, J. A. & López-Jurado, C. Noticiario Ornitológico. Ardeola 66, 205–255 (2019).
Google Scholar 

72.
Logeais, J.-M. Première nidification de l’Élanion blanc Elanus caeruleus en Maine-et-Loire. Crex 13, 45–50 (2015).
Google Scholar 

73.
Balbontín, J., Negro, J. J., Sarasola, J. H., Ferrero, J. J. & Rivera, D. Land-use changes may explain the recent range expansion of the Black-shouldered Kite Elanus caeruleus in southern Europe. Ibis (Lond. 1859). 150, 707–716 (2008).
Google Scholar 

74.
Ferrero, J. J. & Onrubia, A. Elanio azul Elanus caeruleus. In Libro Rojo de las Aves de España (eds Madroño, A. et al.) 113–116 (SEO/BirdLife, Madrid, 2004).
Google Scholar 

75.
Knaus, P. et al. Swiss Breeding Bird Atlas 2013–2016. Distribution and Population Trends of Birds in Switzerland and Liechtenstein (Swiss Ornithological Institute, Sempach, 2018).
Google Scholar 

76.
Pulido-Pastor, A., Márquez, A. L., García-Barros, E. & Real, R. Identification of potential source and sink areas for butterflies on the Iberian Peninsula. Insect Conserv. Divers. 11, 479–492 (2018).
Google Scholar 

77.
Muñoz, A.-R. & Real, R. Assessing the potential range expansion of the exotic monk parakeet in Spain. Divers. Distrib. 12, 656–665 (2006).
Google Scholar 

78.
Niamir, A., Skidmore, A. K., Muñoz, A. R., Toxopeus, A. G. & Real, R. Incorporating knowledge uncertainty into species distribution modelling. Biodivers. Conserv. 28, 571–588 (2019).
Google Scholar 

79.
Muñoz, A.-R. & Real, R. Distribution of Bonelli’s Eagle Aquila fasciata in southern Spain: Scale may matter. Acta Ornithol. 48, 93–101 (2013).
Google Scholar 

80.
Barbosa, A. M. & Real, R. Applying fuzzy logic to comparative distribution modelling: A case study with two sympatric amphibians. Sci. World J. 2012, 1–10 (2012).
Google Scholar 

81.
Maguire, B. J. Niche response structure and the analytical potentials of its relationship to the habitat. Am. Nat. 107, 213–246 (1973).
Google Scholar 

82.
Mill, H. The International Geography (D. Appleton and Company, Boston, 1902).
Google Scholar 

83.
Udvardy, M. D. F. A classification of the Biogeographical Provinces of the World. (1975).

84.
Font, I. Climatología de España y Portugal (Ediciones Universidad de Salamanca, Salamanca, 2000).
Google Scholar 

85.
Crovello, T. J. Quantitative Biogeography: An Overview. Taxon 30, 563–575 (1981).
Google Scholar 

86.
Márquez, A. L., Real, R. & Vargas, J. M. Methods for comparison of biotic regionalizations: The case of pteridophytes in the Iberian Peninsula. Ecography (Cop.) 24, 659–670 (2001).
Google Scholar 

87.
Barreau, D., Bergier, P. & Lesne, L. Lávifaune de l’Oukaïmeden, 2200–3600m (Haut Atlas, Maroc). L’Oiseau la R.F.O.57, 307–367 (1987).

88.
Rodriguez, G., Elorriaga, J. & Ramirez, J. Identification of Atlas Long-Legged Buzzard (Buteo rufinus cirtensis) and its status in Europe. Bird. World 26, 147–173 (2013).
Google Scholar 

89.
IUCN. Buteo rufinus. The IUCN Red List of Threatened Species. Version 2018-2. https://www.iucnredlist.org (2017).

90.
Hijmans, R. J., Cameron, S. E., Parra, J. L., Jones, G. & Jarvis, A. Very high resolution interpolated climate surfaces for global land areas. Int. J. Climatol. 25, 1965–1978 (2005).
Google Scholar 

91.
Kirkevåg, A. et al. Aerosol-cloud-climate interactions in the climate model CAM-Oslo. Tellus Ser. A Dyn. Meteorol. Oceanogr. 60A, 492–512 (2008).
Google Scholar 

92.
Collins, W. J. et al. Development and evaluation of an Earth-System model—HadGEM2. Geosci. Model Dev. 4, 1051–1075 (2011).
ADS  Google Scholar 

93.
Fa, J. E. et al. Integrating sustainable hunting in biodiversity protection in central Africa: Hot spots, weak spots, and strong spots. PLoS ONE 9, e112367 (2014).
ADS  PubMed  PubMed Central  Google Scholar 

94.
Benjamini, Y. & Hochberg, Y. Controlling the False Discovery Rate: A practical and powerful approach to Multiple Testing. J. R. Stat. Soc. Ser. B. 57, 289–300 (1995).
MathSciNet  MATH  Google Scholar 

95.
García, L. V. Controlling the false discovery rate in ecological research. Trends Ecol. Evol. 18, 553–554 (2003).
Google Scholar 

96.
Benjamini, Y. & Yekutieli, D. The control of the false discovery rate in multiple testing under depencency. Ann. Stat. 29, 1165–1188 (2001).
MATH  Google Scholar 

97.
Hosmer, D. W. & Lemeshow, S. Applied Logistic Regression (Wiley, Hoboken, 2000).
Google Scholar 

98.
Real, R., Barbosa, A. M. & Vargas, J. M. Obtaining environmental favourability functions from logistic regression. Environ. Ecol. Stat. 13, 237–245 (2006).
MathSciNet  Google Scholar 

99.
Muñoz, A.-R., Jiménez-Valverde, A., Márquez, A. L., Moleón, M. & Real, R. Environmental favourability as a cost-efficient tool to estimate carrying capacity. Divers. Distrib. 21, 1388–1400 (2015).
Google Scholar 

100.
Hosmer, D. W. & Lemeshow, S. Assessing the fit of the model. In Applied Logistic Regression 157–158 (Wiley, Hoboken, 2000).

101.
Chamorro, D., Nieto, I., Real, R. & Muñoz, A. Wintering areas on the move in the face of warmer winters. Ornis Fenn. 96, 1–14 (2019).
Google Scholar 

102.
Wald, A. Tests of statistical hypotheses concerning several parameters when the number of observations is large. Trans. Am. Math. Soc. 54, 426–482 (1943).
MathSciNet  MATH  Google Scholar 

103.
Lobo, J. M., Jiménez-valverde, A. & Real, R. AUC: A misleading measure of the performance of predictive distribution models. Glob. Ecol. Biogeogr. 17, 145–151 (2008).
Google Scholar 

104.
Romero, D., Olivero, J. & Real, R. Comparative assessment of different methods for using land-cover variables for distribution modelling of Salamandra salamandra longirotris. Environ. Conserv. 40, 48–59 (2012).
Google Scholar 

105.
Cohen, J. A coefficient of agreement for nominal scales. Educ. Psychol. Meas. 20, 37–46 (1960).
Google Scholar 

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

107.
Barbosa, A. M., Real, R., Muñoz, A.-R. & Brown, J. A. New measures for assessing model equilibrium and prediction mismatch in species distribution models. Divers. Distrib. 19, 1333–1338 (2013).
Google Scholar 

108.
Romo, H., García-Barros, E., Márquez, A. L., Moreno, J. C. & Real, R. Effects of climate change on the distribution of ecologically interacting species: Butterflies and their main food plants in Spain. Ecography (Cop.) 37, 1063–1072 (2014).
Google Scholar 

109.
Muñoz, A.-R., Márquez, A. L. & Real, R. Updating known distribution models for forecasting climate change impact on endangered species. PLoS ONE 8, e65462 (2013).
ADS  PubMed  PubMed Central  Google Scholar 

110.
Romero, D., Olivero, J. & Real, R. Accounting for uncertainty in assessing the impact of climate change on biodiversity hotspots in Spain. Anim. Biodivers. Conserv. 42, 355–367 (2019).
Google Scholar 

111.
Daget, P. Ordenation des profils ecologiques. Nat. Monspel. Ser. Bot. 26, 109–128 (1977).
Google Scholar 

112.
Gauch, H. G. Multivariate Analysis in Community Ecology. (Cambridge University Press, Cambridge, 1982). https://doi.org/10.1017/CBO9780511623332.

113.
Real, R., Guerrero, J. C., Antúnez, A., Olivero, J. & Vargas, J. M. Geographic responses of amphibian species to environmental gradients in Southern Spain. I. Individualistic patterns. Boletín la Real Soc. Española Hist. Nat. (Sec. Biol.) 96, 251–261 (2001).
Google Scholar 

114.
Antúnez, A. & Mendonza, M. Factores que determinan el área de distribución geográfica de las especies: Conceptos, modelos y métodos de análisis. In Objetivos y Métodos Biogeográficos. Aplicaciones en Herpetología. Monografías de Herpetología (eds. Vargas, J. M., et al.) vol. 2, 51–72 (Asociación Herpetológica Española, 1992).

115.
Kou, X., Li, Q. & Liu, S. Quantifying species’ range shifts in relation to climate change: A case study of Abies spp. in China. PLoS ONE 6, e23115 (2011).
ADS  CAS  PubMed  PubMed Central  Google Scholar 

116.
US Geological Survey. GTOPO30. (Land processes distributed active archive center (LPDAAC), EROS data center, 1996). More

Study region
The Barrow Peninsula (~1800 km2) is situated on the northern limit of the Arctic Coastal Plain (Fig. 1). The mean annual temperature, precipitation, and snowfall are −11.2 °C, 115 mm, and 958 mm, respectively (1981–2010)35 and the maximum thaw depth ranges from 30 to 90 cm36,37. This continuous permafrost region is characterized by meso-scale (tens to hundreds of square kilometers) drained thaw lake basins (DTLBs) and interstitial tundra9,38, which are composed of a mosaic of fine-scale polygonal tundra landforms (tens to hundreds of square meters). Excluding lakes and rivers, the dominant polygonal tundra landforms in this region includes low-center (LC) polygon, flat-center (FC) polygon, high-center (HC) polygon, coalescent LC polygon, drained slopes (DS), nonpatterned DTLB (nDTLB), and thermokarst ponds, which cover an estimated 34, 24, 16, 11, 11, 3, and 1% of the land surface area, respectively9,13. Due to the similarity in morphological and physiological characteristics of coalescent LC polygons and thermokarst ponds, they are rarely differentiated in field observations. Therefore, both these landforms are combined and referred to as Ponds in the proceeding analysis. Though multiple vegetation communities may be found on each tundra landform, communities typically assemble along a soil moisture gradient representative of each landform21. These community–landform associations are identified as follows: dry Salix heath–DS, dry Luzula heath–HC, moist–wet Carex–Oncophorus meadow–FC, moist–wet Carex–Eriophorum meadow–LC, wet Dupontia meadow–nDTLB, and wet Arctophila pond margin–Pond21.
Model parameterization and validation
We synthesized an extensive collection of field data measured on the Barrow Peninsula to parameterize and validate DOS-TEM (Supplementary Table 1 and Fig. 3). The majority of this data was acquired by scientific initiatives: (1) International Biological Research Program during the early 1970s21,38,39,40, (2) Next Generation Ecosystem Experiments between 2010 and 201641,42,43,44,45,46,47, and (3) Carbon in Arctic Reservoirs Vulnerability Experiment (CARVE) during 2011–201548. In addition, we leveraged key ancillary datasets including: soil carbon pedons (i.e., 100 cm soil cores)22,23,24,30,38,39,49,50,51, vegetation carbon and nitrogen21,38,39,40,52, eddy covariance measurements48, and polygonal tundra landform maps9,13.
Fig. 3: Validation of modeled carbon fluxes and carbon pools.

Monthly net ecosystem exchange (NEE) fluxes measured by the CARVE eddy covariance tower (71°19’22.72”N, 156°35’47.74”W, a) during 2011–2015, were compared with NEE fluxes simulated with DOS-TEM (dashed line in b). Negative NEE indicates carbon uptake, while positive NEE indicates loss. Footprint % indicates the accumulated percentage of measured NEE used to compare with modeled NEE, weighted by polygonal landform (DS drained slope, HC high center, FC flat center, LC low center, nDTLB nonpatterned drained thaw lake basins) using the Kormann and Meixner93 flux footprint model (e.g., a). Modeled carbon pools (colored circles; c) were compared to 44 pedons collected (solid gray circles with standard error bars) and validated against 11 independent random subset of soil carbon pedons (open circles) measured on each respective landform across the Barrow Peninsula.

Full size image

Modeled carbon fluxes were compared to net ecosystem exchange (NEE) measurements from the CARVE tower near Utqiaġvik (71°19′22.72′N, 156°35′47.74′W). The tower footprint (~250 m radius) was located in a heterogeneous tundra site composed of all dominant polygonal tundra landforms (exception of Ponds). Although we identified good correspondence with modeled and measured NEE for most of our observations, DOS-TEM underestimated respiratory losses during the zero-curtain seasonal freeze and thaw isothermal period (e.g., September and October)53, resulting in an underestimate of the 1 to 1 line (R2 = 0.46, p 4 km2), while random error became increasingly positive with scale, increasing by 1.4% for every 1 km2 coarsening of spatial scale. Random errors ranged from 3.9 to 22.8% associated with fine to coarse-scale representation of tundra landforms.
Fig. 6: Misrepresentation of polygonal tundra landforms with scale.

Differences in landform distribution across spatial scales are relative to the highest resolution (0.0009 km2; pie chart). Data are representative of mean landform distributions across six subregions on the Barrow Peninsula. Negative and positive values indicate an overestimate and underestimate of polygonal tundra landforms, respectively.

Full size image

Both the bias error and random error were significantly minimized at fine scales (Fig. 5a), as twenty-first century soil carbon was only overestimated by a maximum of 3.7 and ±7.4%, respectively. This is in contrast to coarser spatial scales as bias and random error sharply increased at 8, 16, and 25 km2 by −6.1% (±10.7%), −17.0% (±22.1%), and −12.6% (±35.5%), respectively (Fig. 5a). The increase in spatial scale led to the overestimation in the area of low productivity thermokarst lakes (1.1% for every 1 km2) and underestimated wet productive landforms such as Ponds (−0.5% for every 1 km2) and LC polygons (−0.5% for every 1 km2; Fig. 6). This underestimation of wet landforms was particularly concerning as wet landforms have been regionally identified as those most sensitive to change59,60,61, while representing a significant proportion of the regional carbon cycle9,60,62,63.
Influence of tundra heterogeneity and model spatial scale
To evaluate the causes, consequences, and mitigation strategies for twenty-first century errors of prediction (i.e., bias and random error), we examined the combined influence of both tundra heterogeneity and model spatial scale. Correlation matrices clarified the potential causes of variable prediction errors, while hierarchical cluster analysis implemented using Euclidean distance and McQuitty linkage methods were used for grouping tundra heterogeneity and model spatial scales with similar errors of prediction to identify potential mitigation strategies or recommendations for future modeling applications.
Correlation matrices supported our presumption that an overestimation of lakes and underestimation of productive wet landforms altered the quantification of landscape-level soil carbon stocks, as bias error was strongly negatively correlated with lake cover (r = −0.98) and positively correlated with wet landforms (r = 0.94; Fig. 7). We found an inverse correlation in bias error as the prevalence of lake cover increased with spatial scale at the expense of nearly all other landforms, but in particular the landforms in low abundance such as tundra ponds (Figs. 6 and 7). Similar to the identified influence of spatial scale on random error (Fig. 4), correlations were highly positively related with model spatial scale (r = 0.99; Fig. 7), reinforcing the impact of coarsening model scale on uncertainty propagation.
Fig. 7: Pearson’s correlations of uncertainty metrics (bias and random error) and spatial attributes.

The larger the bubble the greater the p value. Landform categories dry and wet include spatial data from “DS + HC” and “FC + LC + nDTLB+Pond”, respectively. See Supplementary Fig. 3 for correlation bubble plots of all clusters.

Full size image

Overall, bias error was linked with the misrepresentation of tundra landforms as spatial scale increased (Fig. 7 and Supplementary Fig. 4). Therefore, we next elucidated the influence of heterogeneity and scale on random error. Though random error was correlated with spatial scale, we explored the variability across tundra heterogeneity and scale. The lowest and highest random errors occurred at the finest (≤4 km2) and coarsest (≥16 km2) spatial scales, respectively (Fig. 8). Landform clusters include one or more landforms and landform groups needed to represent tundra heterogeneity on the Barrow Peninsula. Random error was constrained to ±4.5% by considering 5 or 6 tundra landform groups at fine scales. However, at coarse scales these heterogeneous groups also showcased the greatest errors (±28.9%) due to the high number of landforms parameterized within increasingly uncertain landform distributions as scale increased (Figs. 5 and 8). The lowest error among clusters was identified in landform cluster 2 (i.e., ±3.4%; dry and wet), likely due to biogeophysical similarities (i.e., soil anaerobicity, soil available nitrogen, productivity gradients) between dry versus wet landforms (Supplementary Table 1) and similar responses to climate change (e.g., Fig. 4a). Interestingly, even at coarse scales the error found in cluster 2 remained lower than all other landform clusters. Although the “tundra-biome” cluster 1 had a relatively low random error across spatial scales (Fig. 8), this result would not be directly transferable to other modeling applications as we leveraged (i) a robust dataset for model parameterization and (ii) high-resolution polygonal tundra landform maps, currently unavailable across the Arctic for initializing and weighting model parameterization data. The importance of our data assimilation and landform weighting protocol was confirmed by testing the performance of a single unweighted landform parameterization (i.e., HC polygon) extrapolated across the Barrow Peninsula. We found random error to double (±15%) that of cluster 1 at fine scales (≤4 km2) and nearly triple (±45%) at coarse scales ( >8 km2). Therefore, to best simulate dynamically changing carbon pools in permafrost soils, our analysis recommends a minimum of two landform groups (i.e., dry and wet) at a maximum model spatial scale of ≤4 km2 (Fig. 8).
Fig. 8: Heat-map of random error for all tundra heterogeneity and model spatial scales.

Warm to cool colors represent high to low random error (transformed to improve visualization). Hierarchical clustering grouped random error for all landform clusters (i.e., landforms and groups to represent the heterogeneity on the Barrow Peninsula) using a ~50% similarity cut-off for group membership. Mean random errors (transparent white circles) are presented for each landform cluster and model spatial scale.

Full size image

Implications for modeling soil carbon dynamics in Arctic tundra
Current uncertainties among Pan-Arctic model projections reflect inadequate spatial and temporal data needed to initialize, parameterize, and validate key Arctic ecosystem processes55,56,64. This study overcame many of these limitations by leveraging a legacy of data (1973–2016) collected from the data-rich Barrow Peninsula to constrain parameter, climate, and model uncertainties, to improve the representation of Arctic tundra heterogeneity across model spatial scales. We identify a scale-dependent balance between tundra heterogeneity and model spatial scale, linked with the decoupling of actual and simulated tundra landform distributions as spatial scales increased (Figs. 5 and 6). The scale-dependency of model process representation is supported by ground-based assessments, as the drivers of carbon dynamics vary across local (e.g., drainage conditions affecting aerobic/anaerobic processes), regional (e.g., vegetation distribution), and landscape scales (e.g., climate variability). Though we identified relatively minimal differences in carbon accumulation rates between polygonal tundra landforms, this was not necessarily surprising as Arctic coastal tundra landforms are relatively young ( More

1.
Bathurst, R. G. C. Carbonate sediments and their diagenesis (Elsevier, Amsterdam, 1971).
Google Scholar 
2.
Beach, K. et al. Variability in the ecophysiology of Halimeda spp. (Chlorphyta, Bryopsidales) on Conch Reef, Florida Keys, USA. J. Phycol. 39, 633–643 (2003).
CAS  Article  Google Scholar 

3.
Vroom, P. S. et al. Field biology of Halimeda tuna (Bryopsidales, Chlorophyta) across a depth gradient: Comparative growth, survivorship, recruitment, and reproduction. Hydrobiologia 501, 149–166 (2003).
Article  Google Scholar 

4.
Littler, M. M., Littler, D. S., Blair, S. M. & Norris, J. N. Deepest known plant life discovered on an uncharted seamount. Science 227, 57–59 (1985).
ADS  CAS  PubMed  Article  Google Scholar 

5.
Littler, M. M., Littler, D. S., Blair, S. M. & Norris, J. N. Deep-water plant communities from an uncharted seamount off San Salvador Island, Bahamas: Distribution, abundance, and primary productivity. Deep Sea Res. Part A Oceanogr. Res. Pap. 33, 881–892 (1986).
ADS  CAS  Article  Google Scholar 

6.
Blair, S. M. & Norris, J. N. The deep-water species of Halimeda Lamouroux (Halimedaceae, Chlorophyta) from San Salvador Island, Bahamas: Species composition, distribution and depth records. Coral Reefs 6, 227–236 (1988).
ADS  Article  Google Scholar 

7.
Drew, E. A. & Abel, K. M. Studies on Halimeda: II. Reproduction, particularly the seasonality of gametangia formation, in a number of species from the Great Barrier Reef Province. Coral Reefs 6, 207–218 (1988).
ADS  Article  Google Scholar 

8.
Johns, H. D. & Moore, C. H. Reef to basin sediment transport using Halimeda as a sediment tracer, Grand Cayman Island. West Indies. Coral Reefs 6, 187–193 (1988).
ADS  Article  Google Scholar 

9.
Schupp, P. J. & Paul, V. J. Calcium carbonate and secondary metabolites in tropical seaweeds: Variable effects on herbivorous fishes. Ecology 75, 1172–1185 (1993).
Article  Google Scholar 

10.
Littler, M. M. & Littler, D. S. Blade abandonment/proliferation: a novel mechanism for rapid epiphyte control in marine macrophytes. Ecology 80, 1736–1746 (1999).
Article  Google Scholar 

11.
Wefer, G. Carbonate production by algae Halimeda, Penicillus and Padina. Nature 285, 323–324 (1980).
ADS  CAS  Article  Google Scholar 

12.
Wiman, S. K. & McKendree, W. G. Distribution of Halimeda plants and sediments on and around a patch reef near Old Rhodes Key Florida. J. Sediment. Res. 45, 415–421 (1975).
Google Scholar 

13.
Mankiewicz, C. Occurrence and paleocologic significance of Halimeda in late Miocene reefs, southeastern Spain. Coral Reefs 6, 271–279 (1988).
ADS  Article  Google Scholar 

14.
Hillis, L. W. The calcareous reef alga Halimeda (Chlorophyta, Byropsidales): A cretaceous genus that diversified in the cenozoic. Palaeogeogr. Palaeoclimatol. Palaeoecol. 166, 89–100 (2001).
Article  Google Scholar 

15.
Stanley, S. M. Effects of global seawater chemistry on biomineralization: Past, present, and future. Chem. Rev. 108, 4483–4498 (2008).
CAS  PubMed  Article  Google Scholar 

16.
Wefer, G. & Berger, W. H. Stable isotope composition of benthic calcareous algae from Bermuda. J. Sediment. Res. 51, 459–465 (1981).
CAS  Google Scholar 

17.
Flügel, E. Halimeda: paleontological record and palaeoenvironmental significance. Coral Reefs 6, 123–130 (1988).
ADS  Article  Google Scholar 

18.
Rao, V. P. et al. Late quaternary Halimeda bioherms and aragonitic faecal pellet-dominated sediments on the carbonate platform of the western continental shelf of India. Mar. Geol. 121, 293–315 (1994).
ADS  Article  Google Scholar 

19.
Stanley, S. M., Ries, J. B. & Hardie, L. A. Nonlinear partial differential equations and applications: From the Cover: Low-magnesium calcite produced by coralline algae in seawater of Late Cretaceous composition. Proc. Natl. Acad. Sci. 99, 15323–15326 (2002).
ADS  CAS  PubMed  Article  Google Scholar 

20.
Stanley, S. M. Influence of seawater chemistry on biomineralization throughout phanerozoic time: Paleontological and experimental evidence. Palaeogeogr. Palaeoclimatol. Palaeoecol. 232, 214–236 (2006).
Article  Google Scholar 

21.
Brachert, T. C. et al. High salinity variability during the early Messinian revealed by stable isotope signatures from vermetid and Halimeda reefs of the Mediterranean region. Geol. Romana 40, 51–66 (2007).
Google Scholar 

22.
Urey, H. C. The thermodynamic properties of isotopic substances. J. Chem. Soc. 31, 562–581 (1947).
Article  Google Scholar 

23.
Urey, H. C., Lowenstam, H. A., Epstein, S. & McKINNEY, C. R. Measurement of paleotemperatures and temperatures of the Upper Cretaceous of England, Denmark, and the southeastern United States. Geol. Soc. Am. Bull. 62, 399 (1951).
ADS  CAS  Article  Google Scholar 

24.
Epstein, S., Buchsbaum, R., Lowenstam, H. & Urey, H. C. Carbonate-water isotopic temperature scale. Geol. Soc. Am. Bull. 62, 417 (1951).
ADS  CAS  Article  Google Scholar 

25.
Epstein, S., Buchsbaum, R., Lowenstam, H. A. & Urey, H. C. Revised carbonate-water isotopic temperature scale. Geol. Soc. Am. Bull. 64, 1315 (1953).
ADS  CAS  Article  Google Scholar 

26.
Tabouret, H. et al. Simultaneous use of strontium: calcium and barium: calcium ratios in otoliths as markers of habitat: Application to the European eel (Anguilla anguilla) in the Adour basin South West France. Mar. Environ. Res. 70, 35–45 (2010).
CAS  PubMed  Article  Google Scholar 

27.
Ren, L. et al. Deconvolving the δ18O seawater component from subseasonal coral δ18O and Sr/Ca at Rarotonga in the southwestern subtropical Pacific for the period 1726 to 1997. Geochim. Cosmochim. Acta 67, 1609–1621 (2003).
ADS  CAS  Article  Google Scholar 

28.
Grossman, E. L. & Ku, T.-L. Oxygen and carbon isotope fractionation in biogenic aragonite: Temperature effects. Chem. Geol. 59, 59–74 (1986).
CAS  Article  Google Scholar 

29.
Weber, J. N. & Woodhead, P. M. J. Temperature dependence of oxygen-18 concentration in reef coral carbonates. J. Geophys. Res. 77, 463–473 (1972).
ADS  CAS  Article  Google Scholar 

30.
Patterson, W. P. & Walter, L. M. Syndepositional diagenesis of modern platform carbonates: Evidence from isotopic and minor element data. Geology 22, 127–130 (1994).
ADS  CAS  Article  Google Scholar 

31.
Hillis-Colinvaux, L. Ecology and taxonomy of Halimeda: primary producer of coral reefs. In Advances in Marine Biology 1–327 (Academic Press, New york, 1980).
Google Scholar 

32.
Lee, D. & Carpenter, S. J. Isotopic disequilibrium in marine calcareous algae. Chem. Geol. 172, 307–329 (2001).
ADS  CAS  Article  Google Scholar 

33.
Wizemann, A., Meyer, F. W. & Westphal, H. A new model for the calcification of the green macro-alga Halimeda opuntia (Lamouroux). Coral Reefs 33, 951–964 (2014).
ADS  Article  Google Scholar 

34.
Holmes, C. W. delta 18 O variations in the Halimeda of Virgin Islands sands; evidence of cool water in the Northeast Caribbean, late Holocene. J. Sediment. Res. 53, 429–438 (1983).
CAS  Google Scholar 

35.
Peach, K. E., Koch, M. S., Blackwelder, P. L., Guerrero-Given, D. & Kamasawa, N. Primary utricle structure of six Halimeda species and potential relevance for ocean acidification tolerance. Bot. Mar. 60, 1–11 (2017).
Article  Google Scholar 

36.
Peach, K. E., Koch, M. S., Blackwelder, P. L. & Manfrino, C. Calcification and photophysiology responses to elevated pCO2 in six Halimeda species from contrasting irradiance environments on Little Cayman Island reefs. J. Exp. Mar. Biol. Ecol. 486, 114–126 (2017).
CAS  Article  Google Scholar 

37.
Aharon, P. Recorders of reef environment histories: Stable isotopes in corals, giant clams, and calcareous algae. Coral Reefs 10, 71–90 (1991).
ADS  Article  Google Scholar 

38.
Wefer, G. & Berger, W. H. Isotope paleontology: growth and composition of extant calcareous species. Mar. Geol. 100, 207–248 (1991).
ADS  CAS  Article  Google Scholar 

39.
Wefer, G. & Killingley, J. S. Carbon isotopes in organic matter from a benthic alga Halimeda incrassata (Bermuda): Effects of light intensity. Chem. Geol. 59, 321–326 (1986).
CAS  Article  Google Scholar 

40.
McConnaughey, T. 13C and 18O isotopic disequilibrium in biological carbonates: I. Patterns. Geochim. Cosmochim. Acta 53, 151–162 (1989).
ADS  CAS  Article  Google Scholar 

41.
Mcconnaughey, T. 13C and 18O isotopic disequilibrium in biological carbonates: II. In vitro simulation of kinetic isotope effects. Geochim. Cosmochim. Acta 53, 163–171 (1989).
ADS  CAS  Article  Google Scholar 

42.
Robbins, L. L., Knorr, P. O. & Hallock, P. Response of Halimeda to ocean acidification: field and laboratory evidence. Biogeosci. Discuss. 6, 4895–4918 (2009).
ADS  Article  Google Scholar 

43.
Price, N., Hamilton, S., Tootell, J. & Smith, J. Species-specific consequences of ocean acidification for the calcareous tropical green algae Halimeda. Mar. Ecol. Prog. Ser. 440, 67–78 (2011).
ADS  CAS  Article  Google Scholar 

44.
Sinutok, S., Hill, R., Doblin, M. A., Wuhrer, R. & Ralph, P. J. Warmer more acidic conditions cause decreased productivity and calcification in subtropical coral reef sediment-dwelling calcifiers. Limnol. Oceanogr. 56, 1200–1212 (2011).
ADS  CAS  Article  Google Scholar 

45.
Barry, S. C., Frazer, T. K. & Jacoby, C. A. Production and carbonate dynamics of Halimeda incrassata (Ellis) Lamouroux altered by Thalassia testudinum Banks and Soland ex König. J. Exp. Mar. Biol. Ecol. 444, 73–80 (2013).
CAS  Article  Google Scholar 

46.
Vogel, N. et al. Calcareous green alga Halimeda tolerates ocean acidification conditions at tropical carbon dioxide seeps: Halimeda growing at CO 2 seeps. Limnol. Oceanogr. 60, 263–275 (2015).
ADS  CAS  Article  Google Scholar 

47.
Campbell, J. E., Fisch, J., Langdon, C. & Paul, V. J. Increased temperature mitigates the effects of ocean acidification in calcified green algae (Halimeda spp). Coral Reefs 35, 357–368 (2016).
ADS  Article  Google Scholar 

48.
Peach, K., Koch, M. & Blackwelder, P. Effects of elevated pCO2 and irradiance on growth, photosynthesis and calcification in Halimeda discoidea. Mar. Ecol. Prog. Ser. 544, 143–158 (2016).
ADS  CAS  Article  Google Scholar 

49.
Prathep, A., Kaewsrikhaw, R., Mayakun, J. & Darakrai, A. The effects of light intensity and temperature on the calcification rate of Halimeda macroloba. J. Appl. Phycol. 30, 3405–3412 (2018).
CAS  Article  Google Scholar 

50.
Teichberg, M., Fricke, A. & Bischof, K. Increased physiological performance of the calcifying green macroalga Halimeda opuntia in response to experimental nutrient enrichment on a Caribbean coral reef. Aquat. Bot. 104, 25–33 (2013).
CAS  Article  Google Scholar 

51.
Duarte, C. M. et al. Stable Isotope (δ13C, δ15N, δ18O, δD) Composition and Nutrient Concentration of Red Sea Primary Producers. Front. Mar. Sci. 5, 298 (2018).
Article  Google Scholar 

52.
Leichter, J. J., Deane, G. B. & Stokes, M. D. Spatial and temporal variability of internal wave forcing on a coral reef. J. Phys. Oceanogr. 35, 1945–1962 (2005).
ADS  Article  Google Scholar 

53.
Leichter, J. J., Stewart, H. L. & Miller, S. L. Episodic nutrient transport to Florida coral reefs. Limnol. Oceanogr. 48, 1394–1407 (2003).
ADS  Article  Google Scholar 

54.
Leichter, J. J. et al. Nitrogen and oxygen isotopic signatures of subsurface nitrate seaward of the Florida Keys reef tract. Limnol. Oceanogr. 52, 1258–1267 (2007).
ADS  CAS  Article  Google Scholar 

55.
Leichter, J. J., Stokes, M. D., Vilchis, L. I. & Fiechter, J. Regional synchrony of temperature variation and internal wave forcing along the Florida Keys reef tract. J. Geophys. Res. Oceans 119, 548–558 (2014).
ADS  Article  Google Scholar 

56.
Davis, K. A., Leichter, J. J., Hench, J. L. & Monismith, S. G. Effects of western boundary current dynamics on the internal wave field of the Southeast Florida shelf. J. Geophys. Res. 113, C09010 (2008).
ADS  Google Scholar 

57.
Leichter, J. J., Wing, S. R., Miller, S. L. & Denny, M. W. Pulsed delivery of subthermocline water to Conch Reef (Florida Keys) by internal tidal bores. Limnol. Oceanogr. 41, 1490–1501 (1996).
ADS  CAS  Article  Google Scholar 

58.
Lee, T. N. et al. Influence of Florida current, gyres and wind-driven circulation on transport of larvae and recruitment in the Florida Keys coral reefs. Cont. Shelf Res. 12, 971–1002 (1992).
ADS  Article  Google Scholar 

59.
Lee, T. N., Schott, F. A. & Zantopp, R. Florida current: Low-frequency variability as observed with moored current meters during April 1982 to June 1983. Science 227, 298–302 (1985).
ADS  CAS  PubMed  Article  Google Scholar 

60.
Lee, T. N. & Mayer, D. A. Low-frequency current variability and spin-off eddies along the shelf off southeast Florida. J. Mar. Res. 35, 193–220 (1987).
Google Scholar 

61.
Boyer, J. N. & Jones, R. D. A view from the bridge: external and internal forces affecting the ambient water quality of the Florida Keys National Marine Sanctuary (FKNMS). In The Everglades, Florida Bay, and Coral Reefs of the Florida Keys: An Ecosystem Sourcebook (eds Porter, J. W. & Porter, K. G.) 609–628 (CRC Press, Boca Raton, 2002).
Google Scholar 

62.
Kruczynski, W. L. & McManus, F. Water quality concerns in the Florida Keys: Sources effects, and solutions. In The Everglades, Florida Bay, and Coral Reefs of the Florida Keys: An Ecosystem Sourcebook (eds Porter, J. W. & Porter, K. G.) 827–881 (CRC Press, Boca Raton, 2002).
Google Scholar 

63.
Schmidt, G.A. Global seawater oxygen-18 database. NASA GISShttp://data.giss.nasa.gov/o18data/ (1999).

64.
Grossman, E. L. & Ku, T. L. Aragonite-water isotopic paleotemperature scale based on the benthic foraminifer Hoeglundina elegans. Geol. Soc. Am. Abstr. Prog. 13, 464 (1981).
Google Scholar 

65.
Rose, K. A., Roth, B. M. & Smith, E. P. Skill assessment of spatial maps for oceanographic modeling. J. Mar. Syst. 76, 34–48 (2009).
Article  Google Scholar 

66.
Juillet-Leclerc, A. & Schmidt, G. A calibration of the oxygen isotope paleothermometer of coral aragonite from porites. Geophys. Res. Lett. 28, 4135–4138 (2001).
ADS  CAS  Article  Google Scholar 

67.
Böhm, F. et al. Oxygen isotope fractionation in marine aragonite of coralline sponges. Geochim. Cosmochim. Acta 64, 1695–1703 (2000).
ADS  Article  Google Scholar 

68.
Brown, B. E. Coral bleaching: causes and consequences. Coral Reefs 16, S129–S138 (1997).
Article  Google Scholar 

69.
Knowlton, N. & Jackson, J. B. C. The ecology of coral reefs. In Marine Community Ecology (eds Bertness, M. D. et al.) 395–422 (Sinauer Associates Inc., Sunderland, 2001).
Google Scholar 

70.
Wolanski, E. & Pickard, G. Upwelling by internal tides and kelvin waves at the continental shelf break on the Great Barrier Reef. Mar. Freshw. Res. 34, 65 (1983).
Article  Google Scholar 

71.
Wolanski, E. & Hamner, W. M. Topographically controlled fronts in the ocean and their biological influence. Science 241, 177–181 (1988).
ADS  CAS  PubMed  Article  Google Scholar 

72.
Wolanski, E. & Delesalle, B. Upwelling by internal waves, Tahiti, French Polynesia. Cont. Shelf Res. 15, 357–368 (1995).
ADS  Article  Google Scholar 

73.
Wolanski, E. & Deleersnijder, E. Island-generated internal waves at Scott Reef Western Australia. Continent. Shelf Res. 18, 1649–1666 (1998).
ADS  Article  Google Scholar 

74.
Sandstrom, H. & Elliott, J. A. Internal tide and solitons on the Scotian Shelf: A nutrient pump at work. J. Geophys. Res. 89, 6415 (1984).
ADS  Article  Google Scholar 

75.
Holloway, P. E. Internal hydraulic jumps and solitons at a shelf break region on the Australian North West Shelf. J. Geophys. Res. 92, 5405 (1987).
ADS  Article  Google Scholar 

76.
Pineda, J. Predictable upwelling and the shoreward transport of planktonic larvae by internal tidal bores. Science 253, 548–549 (1991).
ADS  CAS  PubMed  Article  Google Scholar 

77.
MacKinnon, J. A. & Gregg, M. C. Mixing on the late-summer New England shelf—Solibores, shear, and stratification. J. Phys. Oceanogr. 33, 1476–1492 (2003).
ADS  Article  Google Scholar 

78.
Smith, J. E., Smith, C. M., Vroom, P. S., Beach, K. L. & Miller, S. Nutrient and growth dynamics of Halimeda tuna on Conch Reef, Florida Keys: Possible influence of internal tides on nutrient status and physiology. Limnol. Oceanogr. 49, 1923–1936 (2004).
ADS  Article  Google Scholar 

79.
Bemis, B. E., Spero, H. J., Bijma, J. & Lea, D. W. Reevaluation of the oxygen isotopic composition of planktonic foraminifera: Experimental results and revised paleotemperature equations. Paleoceanography 13, 150–160 (1998).
ADS  Article  Google Scholar 

80.
Bemis, B. E., Spero, H. J. & Thunell, R. C. Using species-specific paleotemperature equations with foraminifera: A case study in the Southern California Bight. Mar. Micropaleontol. 46, 405–430 (2002).
ADS  Article  Google Scholar 

81.
Wellington, G. M., Dunbar, R. B. & Merlen, G. Calibration of stable oxygen isotope signatures in Galápagos corals. Paleoceanography 11, 467–480 (1996).
ADS  Article  Google Scholar 

82.
Lloyd, R. M. Variations in the oxygen and carbon isotope ratios of florida bay mollusks and their environmental significance. J. Geol. 72, 84–111 (1964).
ADS  CAS  Article  Google Scholar 

83.
Swart, P. K. et al. The stable oxygen and carbon isotopic record from a coral growing in Florida Bay: A 160 year record of climatic and anthropogenic influence. Palaeogeogr. Palaeoclimatol. Palaeoecol. 123, 219–237 (1996).
Article  Google Scholar 

84.
Corbett, D. R., Dillon, K., Burnett, W. & Chanton, J. Estimating the groundwater contribution into Florida Bay via natural tracers, 222 Rn and CH 4. Limnol. Oceanogr 45, 1546–1557 (2000).
ADS  CAS  Article  Google Scholar 

85.
Dillon, K. et al. Groundwater flow and phosphate dynamics surrounding a high discharge wastewater disposal well in the Florida Keys. J. Hydrol. 284, 193–210 (2003).
ADS  CAS  Article  Google Scholar 

86.
Paytan, A. et al. Submarine groundwater discharge: An important source of new inorganic nitrogen to coral reef ecosystems. Limnol. Oceanogr. 51, 343–348 (2006).
ADS  CAS  Article  Google Scholar 

87.
Gonzalez, L. A. & Lohmann, K. C. Carbon and oxygen isotopic composition of Holocene reefal carbonates. Geology 13, 811–841 (1985).
ADS  CAS  Article  Google Scholar 

88.
Shinn, E. Spur and groove formation on the Florida reef tract. J. Sediment. Res. 33, 291–303 (1963).
Google Scholar 

89.
Multer, H. G., Gischler, E., Lundberg, J., Simmons, K. R. & Shinn, E. A. Key Largo limestone revisited: Pleistocene shelf-edge facies, Florida Keys, USA. Facies 46, 229–271 (2002).
Article  Google Scholar 

90.
Lidz, B. H., Reich, C. D., Peterson, R. L. & Shinn, E. A. New maps, new information: Coral reefs of the Florida keys. J. Coastal Res. 222, 260–282 (2006).
Article  Google Scholar 

91.
Deane, G. B. & Stokes, M. D. A robust single-cable sensor array for oceanographic use. IEEE J. Oceanic Eng. 27, 760–767 (2002).
ADS  Article  Google Scholar 

92.
Stokes, M. D., Leichter, J. J., Wing, S. & Frew, R. Temperature variability and algal isotopic heterogeneity on a Floridian coral reef. Mar. Ecol. 32, 364–379 (2011).
ADS  CAS  Article  Google Scholar 

93.
Epstein, S. & Mayeda, T. Variation of O18 content of waters from natural sources. Geochim. Cosmochim. Acta 4, 213–224 (1953).
ADS  CAS  Article  Google Scholar 

94.
Monserud, R. A. & Leemans, R. Comparing global vegetation maps with the Kappa statistic. Ecol. Model. 62, 275–293 (1992).
Article  Google Scholar 

95.
Pontius, R. Quantification error versus location error in comparison of categorical maps. Photogr. Eng. Remote Sens. 66, 1011–1016 (2000).
Google Scholar 

96.
Hagen, A. Fuzzy set approach to assessing similarity of categorical maps. Int. J. Geogr. Inf. Sci. 17, 235–249 (2003).
Article  Google Scholar 

97.
Hagen-Zanker, A., Straatman, B. & Uljee, I. Further developments of a fuzzy set map comparison approach. Int. J. Geogr. Inf. Sci. 19, 769–785 (2005).
Article  Google Scholar 

98.
Hagen-Zanker, A. Map comparison methods that simultaneously address overlap and structure. J. Geograph. Syst. 8, 165–185 (2006).
ADS  Article  Google Scholar 

99.
Visser, H. & de Nijs, T. The map comparison kit. Environ. Modell. Softw. 21, 346–358 (2006).
Article  Google Scholar 

100.
Legendre, P. & Fortin, M. J. Spatial pattern and ecological analysis. Vegetatio 80, 107–138 (1989).
Article  Google Scholar 

101.
Fernandez, M. et al. Locality uncertainty and the differential performance of four common niche-based modelling techniques. Biodiv. Inf. 6, 2 (2009).
Article  Google Scholar 

102.
MATLAB. version 9.4.0 (R2018a). Natick, Massachusetts: The MathWorks Inc. (2018).

103.
Fernandez, M. et al. Locality uncertainty and the differential performance of four common niche-based modeling techniques. Biodiv. Inf. 6, 2 (2009).
Article  Google Scholar  More

1.
Duchene, O., Vian, J.-F. & Celette, F. Intercropping with legume for agroecological cropping systems: complementarity and facilitation processes and the importance of soil microorganisms. A review. Agric. Ecosyst. Environ. 240, 148–161. https://doi.org/10.1016/j.agee.2017.02.019 (2017).
Article  Google Scholar 
2.
Tilman, D., Cassman, K. G., Matson, P. A., Naylor, R. & Polasky, S. Agricultural sustainability and intensive production practices. Nature 418, 671. https://doi.org/10.1038/nature01014 (2002).
ADS  CAS  Article  Google Scholar 

3.
Duchene, O. et al. Integrating multipurpose perennial grains crops in Western European farming systems. Agric. Ecosyst. Environ. 284, 106591. https://doi.org/10.1016/j.agee.2019.106591 (2019).
Article  Google Scholar 

4.
Karamesouti, M. et al. Land-use and land degradation processes affecting soil resources: evidence from a traditional Mediterranean cropland (Greece). CATENA 132, 45–55. https://doi.org/10.1016/j.catena.2015.04.010 (2015).
Article  Google Scholar 

5.
Guo, J. H. et al. Significant acidification in major Chinese croplands. Science 327, 1008–1010. https://doi.org/10.1126/science.1182570 (2010).
ADS  CAS  Article  Google Scholar 

6.
Brinkmann, N. et al. Intensive tropical land use massively shifts soil fungal communities. Sci. Rep. 9, 3403. https://doi.org/10.1038/s41598-019-39829-4 (2019).
ADS  CAS  Article  PubMed  PubMed Central  Google Scholar 

7.
Muhammed, S. E. et al. Impact of two centuries of intensive agriculture on soil carbon, nitrogen and phosphorus cycling in the UK. Sci. Total Environ. 634, 1486–1504. https://doi.org/10.1016/j.scitotenv.2018.03.378 (2018).
ADS  CAS  Article  PubMed  PubMed Central  Google Scholar 

8.
Galloway, J. N. et al. Transformation of the nitrogen cycling: recent trends, questions, and potential solutions. Science 320, 889–892. https://doi.org/10.1126/science.1136674 (2008).
ADS  CAS  Article  Google Scholar 

9.
Zhao, Y. et al. Soil acidification in China: is controlling SO2 emissions enough? Environ. Sci. Technol. 43, 8021–8026. https://doi.org/10.1021/es901430n (2009).
ADS  CAS  Article  Google Scholar 

10.
Alonso, F. P., Arias, J. S., Fernández, R. O., Fernández, P. G. & Serrano, R. E. Agronomic implications of the supply of lime and gypsum by-products to palexerults from western Spain. Soil Sci. 171, 65–81. https://doi.org/10.1097/01.ss.0000200557.253069.50 (2006).
ADS  CAS  Article  Google Scholar 

11.
Chiu, Y.-W. et al. Alterations of biochemical indicators in hepatopancreas of the golden apple snail, Pomacea canaliculata, from paddy fields in Taiwan. J. Environ. Biol. 35, 667–673 (2014).
PubMed  PubMed Central  Google Scholar 

12.
Dong, S., Zheng, G., Yu, X. & Fu, C. Biological control of golden apple snail, Pomacea canaliculata by Chinese soft-shelled turtle, Pelodiscus sinensis in the wild rice, Zizania latifolia field. Sci. Agric. 69, 142–146. https://doi.org/10.1590/S0103-90162012000200009 (2012).
Article  Google Scholar 

13.
Guo, J., Zhang, J.-E., Zhao, B., Luo, M. & Zhang, C. The role of spotted green pufferfish Tetraodon nigroviridis in controlling golden apple snail Pomacea canaliculata: an effective biological control approach involving a new agent. Biocontrol Sci. Technol. 26, 1100–1112. https://doi.org/10.1080/09583157.2016.1185511 (2016).
Article  Google Scholar 

14.
Olivier, H. M., Jenkins, J. A., Berhow, M. & Carter, J. A pilot study testing a natural and a synthetic molluscicide for controlling invasive apple snails (Pomacea maculata). Bull. Environ. Contam. Toxicol. 96, 289–294. https://doi.org/10.1007/s00128-015-1709-z (2016).
CAS  Article  PubMed  PubMed Central  Google Scholar 

15.
Litsinger, J. & Estano, D. B. Management of the golden apple snail Pomacea canaliculata (Lamarck) in rice. Crop Prot. 12, 363–370. https://doi.org/10.1016/0261-2194(93)90079-X (1993).
Article  Google Scholar 

16.
Teo, S. S. Evaluation of different duck varieties for the control of the golden apple snail (Pomacea canaliculata) in transplanted and direct seeded rice. Crop Prot. 20, 599–604. https://doi.org/10.1016/S0261-2194(01)00029-1 (2001).
Article  Google Scholar 

17.
Halwart, M. The golden apple snail Pomacea canaliculata in Asian rice farming systems: present impact and future threat. Int. J. Pest Manag. 40, 199–206. https://doi.org/10.1080/09670879409371882 (1994).
Article  Google Scholar 

18.
Serra, A. The use of golden snail Pomacea Sp. as animal feed in the Philippines. Tropicultura 15, 40–43 (1997).
Google Scholar 

19.
Wang, J. et al. Using golden apple snail to mitigate its invasion and improve soil quality: a biocontrol approach. Environ. Sci. Pollut. Res. 27, 14903–14914. https://doi.org/10.1007/s11356-020-07998-9 (2020).
CAS  Article  Google Scholar 

20.
Haynes, R. J. & Zhou, Y.-F. Effect of pH and added slag on the extractability of Si in two Si-deficient sugarcane soils. Chemosphere 193, 431–437. https://doi.org/10.1016/j.chemosphere.2017.10.175 (2018).
ADS  CAS  Article  PubMed  Google Scholar 

21.
Rousk, J. et al. Soil bacterial and fungal communities across a pH gradient in an arable soil. ISME J. 4, 1340–1351. https://doi.org/10.1038/ismej.2010.58 (2010).
Article  PubMed  Google Scholar 

22.
Lauber, C. L., Strickland, M. S., Bradford, M. A. & Fierer, N. The influence of soil properties on the structure of bacterial and fungal communities across land-use types. Soil Biol. Biochem. 40, 2407–2415. https://doi.org/10.1016/j.soilbio.2008.05.021 (2008).
CAS  Article  Google Scholar 

23.
Ramirez, K. S., Craine, J. M. & Fierer, N. Consistent effects of nitrogen amendments on soil microbial communities and processes across biomes. Glob. Change Biol. 18, 1918–1927. https://doi.org/10.1111/j.1365-2486.2012.02639.x (2012).
ADS  Article  Google Scholar 

24.
Campbell, B. J., Polson, S. W., Hanson, T. E., Mack, M. C. & Schuur, E. A. The effect of nutrient deposition on bacterial communities in Arctic tundra soil. Environ. Microbiol. 12, 1842–1854. https://doi.org/10.1111/j.1462-2920.2010.02189.x (2010).
CAS  Article  PubMed  Google Scholar 

25.
Yu, H. et al. Responses of soil biological traits and bacterial communities to nitrogen fertilization mediate maize yields across three soil types. Soil Till. Res. 185, 61–69. https://doi.org/10.1016/j.still.2018.08.017 (2019).
Article  Google Scholar 

26.
Fierer, N. & Jackson, R. B. The diversity and biogeography of soil bacterial communities. Proc. Natl. Acad. Sci. 103, 626–631. https://doi.org/10.1073/pnas.0507535103 (2006).
ADS  CAS  Article  PubMed  Google Scholar 

27.
Roesch, L. F. W. et al. Pyrosequencing enumerates and contrasts soil microbial diversity. ISME J. 1, 283–290. https://doi.org/10.1038/ismej.2007.53 (2007).
CAS  Article  PubMed  PubMed Central  Google Scholar 

28.
Waring, B. G., Weintraub, S. R. & Sinsabaugh, R. L. Ecoenzymatic stoichiometry of microbial nutrient acquisition in tropical soils. Biogeochemistry 117, 101–113. https://doi.org/10.1007/s10533-013-9849-x (2014).
CAS  Article  Google Scholar 

29.
Baldrian, P. Microbial enzyme-catalyzed processes in soils and their analysis. Plant Soil Environ. 55, 370–378 (2009).
CAS  Article  Google Scholar 

30.
Jian, S. et al. Soil extracellular enzyme activities, soil carbon and nitrogen storage under nitrogen fertilization: a meta-analysis. Soil Biol. Biochem. 101, 32–43 (2016).
CAS  Article  Google Scholar 

31.
Frankenberger, W. & Bingham, F. Influence of salinity on soil enzyme activities 1. Soil Sci. Soc. Am. J. 46, 1173–1177. https://doi.org/10.2136/sssaj1982.03615995004600060011x (1982).
ADS  CAS  Article  Google Scholar 

32.
Wang, Y. et al. The influence of soil heavy metals pollution on soil microbial biomass, enzyme activity, and community composition near a copper smelter. Ecotoxicol. Environ. Saf. 67, 75–81. https://doi.org/10.1016/j.ecoenv.2006.03.007 (2007).
CAS  Article  Google Scholar 

33.
Camenzind, T., Hättenschwiler, S., Treseder, K. K., Lehmann, A. & Rillig, M. C. Nutrient limitation of soil microbial processes in tropical forests. Ecol. Monogr. 88, 4–21. https://doi.org/10.1002/ecm.1279 (2018).
Article  Google Scholar 

34.
Alvarez, G. et al. Catalytic power of enzymes decreases with temperature: New insights for understanding soil C cycling and microbial ecology under warming. Glob. Change Biol. 24, 4238–4250. https://doi.org/10.1111/gcb.14281 (2018).
ADS  Article  Google Scholar 

35.
Semer, R. & Reddy, K. R. Evaluation of soil washing process to remove mixed contaminants from a sandy loam. J. Hazard. Mater. 45, 45–57. https://doi.org/10.1016/0304-3894(96)82887-1 (1996).
CAS  Article  Google Scholar 

36.
He, Z., Baligar, V., Ritchey, K., Martens, D. & Kemper, W. Factors affecting phosphate rock dissolution in acid soil amended with liming materials and cellulose. Soil Sci. Soc. Am. J. 60, 1596–1601. https://doi.org/10.2136/sssaj1996.03615995006000050045x (1996).
ADS  CAS  Article  Google Scholar 

37.
Moreno, J., Hernández, T. & Garcia, C. Effects of a cadmium-contaminated sewage sludge compost on dynamics of organic matter and microbial activity in an arid soil. Biol. Fertil. Soils 28, 230–237. https://doi.org/10.1007/s003740050487 (1999).
CAS  Article  Google Scholar 

38.
Lauber, C. L., Zhou, N., Gordon, J. I., Knight, R. & Fierer, N. Effect of storage conditions on the assessment of bacterial community structure in soil and human-associated samples. FEMS Microbiol. Lett. 307, 80–86. https://doi.org/10.1111/j.1574-6968.2010.01965.x (2010).
CAS  Article  PubMed  PubMed Central  Google Scholar 

39.
Dong, L. et al. Enhanced nitrate-nitrogen removal by modified attapulgite-supported nanoscale zero-valent iron treating simulated groundwater. J. Environ. Manag. 213, 151–158. https://doi.org/10.1016/j.jenvman.2018.02.073 (2018).
CAS  Article  Google Scholar 

40.
Bray, R. H. & Kurtz, L. Determination of total, organic, and available forms of phosphorus in soils. Soil Sci. 59, 39–46 (1945).
ADS  CAS  Article  Google Scholar 

41.
Au-Bell, C. W. et al. High-throughput fluorometric measurement of potential soil extracellular enzyme activities. J. Vis. Exp. 81, e50961. https://doi.org/10.3791/50961 (2013).
CAS  Article  Google Scholar 

42.
Edgar, R. C. UPARSE: highly accurate OTU sequences from microbial amplicon reads. Nat. Methods 10, 996. https://doi.org/10.1038/nmeth.2604 (2013).
CAS  Article  Google Scholar 

43.
Edgar, R. C., Haas, B. J., Clemente, J. C., Quince, C. & Knight, R. UCHIME improves sensitivity and speed of chimera detection. Bioinformatics 27, 2194–2200. https://doi.org/10.1093/bioinformatics/btr381 (2011).
CAS  Article  PubMed  PubMed Central  Google Scholar 

44.
Quast, C. et al. The SILVA ribosomal RNA gene database project: improved data processing and web-based tools. Nucleic Acids Res. 41, D590–D596. https://doi.org/10.1093/nar/gks1219 (2012).
CAS  Article  PubMed  PubMed Central  Google Scholar 

45.
Meenatchi, R. et al. Revealing the impact of global mass bleaching on coral microbiome through 16S rRNA gene-based metagenomic analysis. Microbiol. Res. 233, 126408. https://doi.org/10.1016/j.micres.2019.126408 (2020).
CAS  Article  Google Scholar 

46.
Gotelli, N. J. & Colwell, R. K. Estimating species richness. Biol. Divers. Front. Meas. Assess. 12, 39–54 (2011).
Google Scholar 

47.
Keenan, S. W., Schaeffer, S. M., Jin, V. L. & DeBruyn, J. M. Mortality hotspots: nitrogen cycling in forest soils during vertebrate decomposition. Soil Biol. Biochem. 121, 165–176. https://doi.org/10.1016/j.soilbio.2018.03.005 (2018).
CAS  Article  Google Scholar 

48.
Birla, A., Singh, B., Upadhyay, S. & Sharma, Y. Kinetics studies of synthesis of biodiesel from waste frying oil using a heterogeneous catalyst derived from snail shell. Biores. Technol. 106, 95–100. https://doi.org/10.1016/j.biortech.2011.11.065 (2012).
CAS  Article  Google Scholar 

49.
Lauber, C. L., Hamady, M., Knight, R. & Fierer, N. Pyrosequencing-based assessment of soil pH as a predictor of soil bacterial community structure at the continental scale. Appl. Environ. Microbiol. 75, 5111–5120. https://doi.org/10.1128/AEM.00335-09 (2009).
CAS  Article  PubMed  PubMed Central  Google Scholar 

50.
Nicol, G. W., Leininger, S., Schleper, C. & Prosser, J. I. The influence of soil pH on the diversity, abundance and transcriptional activity of ammonia oxidizing archaea and bacteria. Environ. Microbiol. 10, 2966–2978. https://doi.org/10.1111/j.1462-2920.2008.01701.x (2008).
CAS  Article  PubMed  PubMed Central  Google Scholar 

51.
Delgado-Baquerizo, M. et al. A global atlas of the dominant bacteria found in soil. Science 359, 320–325. https://doi.org/10.1126/science.aap9516 (2018).
ADS  CAS  Article  PubMed  PubMed Central  Google Scholar 

52.
Li, Q. et al. Biochar amendment decreases soil microbial biomass and increases bacterial diversity in Moso bamboo (Phyllostachys edulis) plantations under simulated nitrogen deposition. Environ. Res. Lett. 13, 1–10. https://doi.org/10.1088/1748-9326/aab53a (2018).
CAS  Article  Google Scholar 

53.
Xiong, J. et al. Geographic distance and pH drive bacterial distribution in alkaline lake sediments across Tibetan Plateau. Environ. Microbiol. 14, 2457–2466. https://doi.org/10.1111/j.1462-2920.2012.02799.x (2012).
CAS  Article  PubMed  PubMed Central  Google Scholar 

54.
Leff, J. W. et al. Consistent responses of soil microbial communities to elevated nutrient inputs in grasslands across the globe. Proc. Natl. Acad. Sci. 112, 10967–10972. https://doi.org/10.1073/pnas.1508382112 (2015).
ADS  CAS  Article  PubMed  PubMed Central  Google Scholar 

55.
Gomez, E., Ferreras, L. & Toresani, S. Soil bacterial functional diversity as influenced by organic amendment application. Biores. Technol. 97, 1484–1489. https://doi.org/10.1016/j.biortech.2005.06.021 (2006).
CAS  Article  Google Scholar 

56.
Freitag, T. E., Chang, L., Clegg, C. D. & Prosser, J. I. Influence of inorganic nitrogen management regime on the diversity of nitrite-oxidizing bacteria in agricultural grassland soils. Appl. Environ. Microbiol. 71, 8323–8334. https://doi.org/10.1128/AEM.71.12.8323-8334.2005 (2005).
CAS  Article  PubMed  PubMed Central  Google Scholar 

57.
Zhang, Y. et al. Fertilization shapes bacterial community structure by alteration of soil pH. Front. Microbiol. 8, 1325. https://doi.org/10.3389/fmicb.2017.01325 (2017).
Article  PubMed  PubMed Central  Google Scholar 

58.
Nemergut, D. R. et al. The effects of chronic nitrogen fertilization on alpine tundra soil microbial communities: implications for carbon and nitrogen cycling. Environ. Microbiol. 10, 3093–3105. https://doi.org/10.1111/j.1462-2920.2008.01735.x (2008).
CAS  Article  PubMed  PubMed Central  Google Scholar 

59.
Zhang, H. et al. Gemmatimonas aurantiaca gen. nov., sp. nov., a Gram-negative, aerobic, polyphosphate-accumulating micro-organism, the first cultured representative of the new bacterial phylum Gemmatimonadetes phyl. nov. Int. J. Syst. Evol. Microbiol. 53, 1155–1163. https://doi.org/10.1099/ijs.0.02520-0 (2003).
CAS  Article  PubMed  PubMed Central  Google Scholar 

60.
Sheng, Y. & Zhu, L. Biochar alters microbial community and carbon sequestration potential across different soil pH. Sci. Total Environ. 622, 1391–1399. https://doi.org/10.1016/j.scitotenv.2017.11.337 (2018).
ADS  CAS  Article  PubMed  Google Scholar 

61.
Cobo-Díaz, J. F. et al. Metagenomic assessment of the potential microbial nitrogen pathways in the rhizosphere of a Mediterranean forest after a wildfire. Microb. Ecol. 69, 895–904. https://doi.org/10.1007/s00248-015-0586-7 (2015).
Article  PubMed  Google Scholar 

62.
Wang, C., Liu, D. & Bai, E. Decreasing soil microbial diversity is associated with decreasing microbial biomass under nitrogen addition. Soil Biol. Biochem. 120, 126–133. https://doi.org/10.1016/j.soilbio.2018.02.003 (2018).
CAS  Article  Google Scholar 

63.
Fierer, N., Bradford, M. A. & Jackson, R. B. Toward an ecological classification of soil bacteria. Ecology 88, 1354–1364. https://doi.org/10.1890/05-1839 (2007).
Article  PubMed  Google Scholar 

64.
Jones, R. T. et al. A comprehensive survey of soil acidobacterial diversity using pyrosequencing and clone library analyses. ISME J. 3, 442–453. https://doi.org/10.1038/ismej.2008.127 (2009).
CAS  Article  PubMed  PubMed Central  Google Scholar 

65.
Baker, B. J. et al. Community transcriptomic assembly reveals microbes that contribute to deep-sea carbon and nitrogen cycling. ISME J. 7, 1962–1973. https://doi.org/10.1038/ismej.2013.85 (2013).
CAS  Article  PubMed  PubMed Central  Google Scholar 

66.
Xu, N., Tan, G., Wang, H. & Gai, X. Effect of biochar additions to soil on nitrogen leaching, microbial biomass and bacterial community structure. Eur. J. Soil Biol. 74, 1–8. https://doi.org/10.1016/j.ejsobi.2016.02.004 (2016).
CAS  Article  Google Scholar 

67.
Prasanna, S. N. R. Soil pH and its role in cyanobacterial abundance and diversity in rice field soils. Appl. Ecol. Environ. Res. 5, 103–113. https://doi.org/10.15666/aeer/0502_103113 (2007).
Article  Google Scholar 

68.
Shen, C. et al. Soil pH drives the spatial distribution of bacterial communities along elevation on Changbai Mountain. Soil Biol. Biochem. 57, 204–211. https://doi.org/10.1016/j.soilbio.2012.07.013 (2013).
CAS  Article  Google Scholar 

69.
Waldrop, M., Balser, T. & Firestone, M. Linking microbial community composition to function in a tropical soil. Soil Biol. Biochem. 32, 1837–1846. https://doi.org/10.1016/S0038-0717(00)00157-7 (2000).
CAS  Article  Google Scholar 

70.
Bombeo-Tuburan, I., Fukumoto, S. & Rodriguez, E. Use of the golden apple snail, cassava, and maize as feeds for the tiger shrimp, Penaeus monodon, in ponds. Aquaculture 131, 91–100 (1995).
Article  Google Scholar 

71.
DeForest, J. L., Zak, D. R., Pregitzer, K. S. & Burton, A. J. Atmospheric nitrate deposition, microbial community composition, and enzyme activity in northern hardwood forests. Soil Sci. Soc. Am. J. 68, 132–138. https://doi.org/10.2136/sssaj2004.1320 (2004).
ADS  CAS  Article  Google Scholar 

72.
Keeler, B. L., Hobbie, S. E. & Kellogg, L. E. Effects of long-term nitrogen addition on microbial enzyme activity in eight forested and grassland sites: implications for litter and soil organic matter decomposition. Ecosystems 12, 1–15. https://doi.org/10.1007/s10021-008-9199-z (2009).
CAS  Article  Google Scholar 

73.
Carreiro, M., Sinsabaugh, R., Repert, D. & Parkhurst, D. Microbial enzyme shifts explain litter decay responses to simulated nitrogen deposition. Ecology 81, 2359–2365. https://doi.org/10.1890/0012-9658(2000)081[2359:MESELD]2.0.CO;2 (2000).
Article  Google Scholar 

74.
Demisie, W. & Zhang, M. Effect of biochar application on microbial biomass and enzymatic activities in degraded red soil. Afr. J. Agric. Res. 10, 755–766. https://doi.org/10.5897/AJAR2013.8209 (2015).
CAS  Article  Google Scholar 

75.
Zheng, M., Huang, J., Chen, H., Wang, H. & Mo, J. Responses of soil acid phosphatase and beta-glucosidase to nitrogen and phosphorus addition in two subtropical forests in southern China. Eur. J. Soil Biol. 68, 77–84. https://doi.org/10.1016/j.ejsobi.2015.03.010 (2015).
CAS  Article  Google Scholar  More