Ecology

Study site
The study site (N 38° 49′, E 115° 26′) is located in the Guanzhuang Village in Baoding City of Hebei Province, China, in the humid temperate and monsoon climatic zone with the average annual air temperature of 13 °C, annual rainfall of 500 mm, and frost-free period is 210 days. Although the experiment was a one-year, the distribution of precipitation (488.50 mm) and temperature (13.45 °C) during the experimental period (2014–2015) were close to the the latest 10-year averages (2005–2015) (500.19 mm and 13.61 °C) (Fig. 1). Determine the basic nutrient indexes of the 0–20 cm surface soil in the test plot. The soil type is silty loam, consists with 22.55% sand, 71.09% silt and 6.36% clay. Analysis of soil basic characteristics showed that it has a pH of 8.3, and its content of organic matter, total N, available phosphorus (P) and available potassium (K) was 11.27 g kg−1, 1.47 g kg−1, 25.49 mg kg−1 and 127.43 mg kg−1, respectively.
Figure 1

Monthly precipitation and average temperature during the experimental year (2014–2015) and the mean values in the last ten years (2005–2015) in the test area. Data was from meteorological station.

Full size image

Experiment materials
The planting mode of the experimental site was a winter wheat-summer maize rotation, the winter wheat variety was ‘Jinnong 6’ that verage thousand weight was 47.6 g. The summer maize variety was ‘Zhengdan 958′ that average thousand weight 330 g.
Test fertilizers include inorganic fertilizer, organic fertilizer, soil conditioner, compound bateria, amino acid liquid fertilizer and nutrient agent. Inorganic N, P and K in the tested fertilizers were provided by urea (N 46%), superphosphate (P2O5,16%) and potassium chloride (K2O, 54%), respectively, as well as in the form of zinc and humic acid urea which is mainly a combination of N with humic acid (N 46% and HA 1.2%). The organic fertilizer used in the experiment was mainly decomposed chicken manure. But the N content in the chicken manure is 1.32% in wheat season and 4.48% in maize sason. Soil conditioner mainly containing calcium (Ca) and magnesium (Mn), Compound bacteria could fix N potentially, promote root growth, decompose cellulose lignin and thus rapidly to degrade. The number of living bacteria reached 2 billion per gram. Amino acid liquid fertilizer and nutrient agent sprayed according to crop growth to provide the required amino acids and trace elements for plant growth.
Experiment design
Field experiment consisted of five treatments with 3 replicates. The experiment uses a completely randomized block setting, the plot size was 79.2 m2 (13.2 m × 6 m). Before the experiment, no crops were planted in the area and it was idle for more than one year. The five treatments were: CK (zero N), FN (farmers’ traditional inorganic N rate, through mass surveys on actual production), RN (recommended inorganic-N rate, according to the experimental results of many scholars, combined with the local soil N supply, crop straw returning in the previous season as well as wheat or maize N demand for target yield in the current season)22,23, HAN (zinc and humic acid urea, the N supply same as RN), RN40% + HOM (40% inorganic N rate of RN (RN40%) with homemade organic matters (HOM). HOM was an organic control measure, it including organic fertilizer, soil conditioner compound bacteria, amino acid liquid fertilizer and nutrient agent. these constituents and amount according to Shu et al24 (Table 1).
Table 1 Rates of pure N and organic matters in different fertilization treatments.
Full size table

For wheat, N fertilizer was broadcast for ratio of 4:3 (basal to topdressing) in RN40% + HOM, whereas for the rest N treatments the ratio was 1:1. During maize planting, N ratio (basal to topdressing) for all treatments was 2:3.The N, P and K fertilizers for wheat were applied in the form of urea, single superphosphate and potassium chloride, respectively. The amount of N fertilizer applied in different treatments of different crops is different, specific application amount reference Table 1. Except for treatment RN40% + HOM, all treatments have the same amount of single superphosphate and potassium chloride. Single superphosphate (120 kg P2O5 ha−1)and potassium chloride (150 kg K2O ha−1)were used in winter wheat season. For maize, single superphosphate (90 kg P2O5 ha−1) and potassium chloride (150 kg K2O ha−1) were used. P and K fertilizers were applied once before sowing. For the doses of P and K in RN40% + HOM brought by organic fertilizer were firstly assessed (48.4 kg P2O5 ha−1 and 149.3 kg K2O ha−1 for winter wheat; 87.2 P2O5 ha−1 and 19.1 kg K2O ha−1 for maize), remaining amounts were supplemented with chemical P and K fertilizers.
Wheat at a rate of 187.5 kg ha−1 with a row space of 15 cm, was sown on 12 October 2014 and harvested at 7 June 2015. Then, at the same wheat plot, Maize of 37.5 kg ha−1 with a row space of 57 cm, was sown on 18 June 2015 and harvested at 5 October 2015.
N2O sampling and measurements
N2O gas was collected using a closed static chamber from sowing to harvest of wheat and maize25. The sampling box was divided into two parts and made by PVC material: a box body and a base. The upper part of the box body was provided with a gas sampling port sealed with a rubber plug, and a thermometer probe was arranged inside the box body to monitor the soil surface temperature. The box body is 15 cm high and the bottom diameter is 25 cm. The base was ring shaped, and buried into the soil. Gas collection was performed from 9:00 to10:00 am. A 30 mL of air sample was collected at 0, 8, 16 min after closure26. The air samples were taken once at an interval of 7 days in general and subsequently continuous 5 days following fertilization or precipitation. Continue to collect gas samples from the beginning of the experiment. No gas samples are collected during the freezing of wheat field soil in February and March every year. At the same time, the air temperature was measured by a thermometer and the soil moisture in the 0–5 cm depth was measured by soil moisture tester (TK3-BASIC). Gas concentrations were analyzed by using a gas chromatography (Agilent 7890 A, USA), fitted with a 4 mm by 3 m stainless steel column packed with Porapack Q and N2 was used as the carrier gas. The column and the detector temperatures were set at 70 °C and 300 °C, respectively. The standard N2O was supplied from National Center of Standard Measurement.
N2O flux was calculated using the following equation (Wang et al.27).

$${text{F}} = rho times {text{H}} times T_{0} frac{{left( {c_{2} /T_{2} – c_{1} /T_{1} } right)}}{Delta t}$$
(1)

where F is N2O emission flux, ρ = m/22.414, ρ is the density of gas in airtight box, m is molecular weight, H is the height of the static chamber, T0 is 273 K, c1 and c2 are the gas concentration in time of t1and t2, respectively, T1and T2 are gas temperatures, ∆t = t2 − t1 ,where t2 and t1are times.
Cumlative N2O emissions were from the growth season was calculated by the equation:

$${text{T = }}sum {left[ {{{left( {{text{F}}_{{text{i + 1}}} {text{ + F}}_{{text{i}}} } right)} mathord{left/ {vphantom {{left( {{text{F}}_{{text{i + 1}}} {text{ + F}}_{{text{i}}} } right)} {2}}} right. kern-nulldelimiterspace} {2}}} right]} times left( {{text{D}}_{{text{i + 1}}} – {text{D}}_{{text{i}}} } right) times {{{24}} mathord{left/ {vphantom {{{24}} {{1000} times {667} times {15} times {10}^{{ – {6}}} }}} right. kern-nulldelimiterspace} {{1000} times {667} times {15} times {10}^{{ – {6}}} }}$$
(2)

where T is the total amount of N2O emissions from the growth stage (kg N ha−1), Fi and Fi+1 denote the N2O flux of the i and i + 1 sub-sampling (μg N m−2 h−1); Di and Di+1 represent sampling days (d)26.
N2O emission coefficient (EF) was estimated with equation28:

$$EF(% ) = left[ {left( {{text{Cumulative}};{text{N}}_{{2}} {text{O}};{text{emissions}};{text{from}};{text{fertilized}};{text{plots}} – {text{control}};{text{plots}}} right)/{text{N}};{text{fertilizer}};{text{rate}}} right] times 100$$
(3)

Soil sampling and measurements
At wheat and maize maturity, soil samples were collected from depths of 0–20, 20–40 and 40–60 cm with a hand probe from three places in central rows of each plot and mixed together. Fresh soil samples were sieved through a 2 mm, extracted with 1 mol L−1 KCl and a soil-solution ratio of 1:10, and analyzed for inorganic N (mainly including NH4+–N and NO3−–N) contents with continuous flow analysis technique(AA3-HR, Germany)22. Soil moisture and density of each soil layer were measured simultaneously, and the soil residual N in 0–60 cm was calculated. Other soil sample was air-dried and sieved, organic matter and total N content were measured by agrochemical analysis method29.
Plant harvest
For wheat, plants with double rows (1 m length) in each plot were harvested and 20 spikes were selected to count the numbers of effective spikes. All the harvested plant samples were separated into straw (including stem, leaves and remaining of ears) and grains, and the grain yield was calculated to 12.5% moisture content (PM-8188, Japan). Three samples were chosen from each plot and weighted to get the average 1000-grain weight.
For maize, two representative plants in each plot were harvested and separated into straw (including stems, leaves, tassels, husks, cobs) and grains in the central rows. Moreover, 20 ears were continuously selected to thresh and measured grain yield. Grain yield was calculated to 14% moisture content (PM-8188, Japan).
All harvested wheat and maize samples were dried, weighed, ground into powder to measure the total N content using H2SO4–H2O2 Kjeldahl digestion method29.
N balance and N efficiencies
Total N input was comprised of N fertilizer, the initial inorganic N in soil before planting (including both NO3−–N and NH4+–N), pre-crop N straw return (no straw was returned when sowing wheat and the N uptake in maize stage was calculated from pre-wheat straw), N deposition from dry and wet atmosphere and mineralized N in soil. Atmospheric N deposition was derived from Research result by Liu et al.30. N output was comprised of crop uptake, post-harvest residual soil N and apparent N loss. This study calculated soil N to a depth of 0–60 cm. Mineralized soil N, apparent N loss, Nitrogen production effiency (NPE) , Nitrogen agronomic effiency (NAE) and Nitrogen use efficiency (NUE) were calculated as follows:

$$begin{aligned} {text{Mineralized}};{text{N}}left( {{text{kg}},{text{ha}}^{ – 1} } right) & = {text{Crop}};{text{N}};{text{uptake}};{text{in}};{text{CK}} + {text{Post – harvest}};{text{residual}};{text{soil}};{text{N}};{text{in}};{text{CK}}{-}{text{Pre – planting}};{text{soil}};{text{N}};{text{in}};{text{CK}} \ & quad {-}{text{N}};{text{deposition}};{text{from}};{text{atmosphere}};{text{in}};{text{CK}} – {text{Pre – crop}};{text{straw}};{text{return}};{text{N}};{text{in}};{text{CK}} \ end{aligned}$$
(4)

$${text{Apparent}};{text{N}};{text{loss}}left( {{text{kg}},{text{ha}}^{ – 1} } right) = {text{Total}};{text{N}};{text{input}} – {text{crop}};{text{N}};{text{Uptake}} – {text{post – harvest}};{text{residual}};{text{soil}};{text{N}}$$
(5)

$${text{NPE}}left( {{text{kg}},{text{kg}}^{ – 1} } right) = {text{Plant}};{text{yield}}/{text{N}};{text{fertlizer}};{text{rate}}$$
(6)

$${text{NAE}}left( {{text{kg}},{text{kg}}^{ – 1} } right) = left[ {left( {{text{Plant}};{text{yield}};{text{with}};{text{N}};{text{application}} – {text{plant}};{text{yield}};{text{without}};{text{N}};{text{fertlizer}}} right)/{text{N}};{text{fertlizer}};{text{rate}}} right]$$
(7)

$${text{NUE}}left( % right) = left[ {left( {{text{N}};{text{content}};{text{in}};{text{plant}};{text{with}};{text{N}};{text{fertilizer}}{-}{text{N}};{text{content}};{text{in}};{text{plant}};{text{without}};{text{N}};{text{fertilizer}}} right)/{text{N}};{text{fertlizer}};{text{rate}}} right] times 100$$
(8)

Net income analyses
Prices of fertilizers and grains as well as other costs in Chinese Yuan (RMB: 1 USD = 6.71 RMB in the experiment year) were based on local prices. Net income was calculated by the equation:

$${text{Net}};{text{income}} = {text{Output}};{text{value}}{-}{text{fertilizer}};{text{cost}}{-}{text{other}};{text{field}};{text{management}};{text{costs}}$$
(9)

$${text{Output}};{text{value}} = {text{Grain}};{text{yield}} times {text{grain}};{text{price}}$$
(10)

where Fertilizer costs were composed of the prices of inorganic N (3.9 RMB kg−1), P2O5 (5.65 RMB kg−1), K2O (6.5 RMB kg−1), pure N in zinc and humic acid urea (5.0 RMB kg−1), decomposed chicken manure (0.5 RMB kg−1), soil conditioner (2.8 RMB kg−1), compound bacteria, amino acid liquid fertilizer and nutrient agent together (30 RMB kg−1). Other field management costs included seed, labor for fertilization, irrigation, mechanical sowing, etc. Grain prices of wheat and maize during the experiment were 2.2 and 1.8 RMB kg−1, respectively.
Statistical analysis
This research adopted SPSS Statistics 20.0 software (SPSS Inc., Chicago, IL, USA) to date analysis. Through least significant differences (LSD) method, the statistically significant differences were calculated. The differences level was prominent when P  More

Jaw shape variation and diet in small mammals
Using 2D geometric morphometrics (Fig. 2a), we found that jaw shape is a good proxy for diet among small extant mammals. In Fig. 3, taxa with negative PC1 scores have shorter jaws, and taxa with positive PC1 scores have longer jaws; taxa with positive PC2 scores have taller ascending rami and taxa with negative PC2 scores have shorter ascending rami. Among extant mammals, most dietary categories (excluding omnivores) can be distinguished along PC1 (Fig. 3a): herbivores plot at the negative end of PC1, insectivores towards the positive end, and carnivores in between. These categories are also statistically different from each other (Table 2), showing that jaw shape can distinguish between most major dietary types. However, our data cannot distinguish between carnivores and omnivores.
Fig. 2: Data acquired from the jaws of Mesozoic and extant small mammals.

a Jaw landmarking regimen used in this study. Modified from ref. 12. In orange: six fixed landmarks; in blue: 58 sliding semi landmarks. b Moment arm measurements taken in this study. Modified from ref. 19.

Full size image

Fig. 3: Scatter plots of the principal component analysis (PCA) results (PC1 vs. PC2).

a Extant taxa, b extinct taxa. Convex hulls shown for extant insectivores (yellow), carnivores (red), omnivores (purple) and herbivores (blue). Icon colors indicate known dietary categories of extant mammals and suggested dietary categories for Mesozoic mammals (obtained from the literature). See Table 1 for taxon names.

Full size image

Table 2 Summary of the Procrustes ANOVA (Type II, Conditional SS) performed for jaw shape data as a function of dietary group.
Full size table

Data on the jaw shape of Mesozoic mammals were projected onto the extant taxa morphospace (Fig. 3b). In order to determine whether jaw shape could be used as a dietary proxy in Mesozoic mammals, we obtained previous independent determinations of likely diets, which variously employed dental morphology, tooth wear facets and body size (e.g., see refs. 1,7,12,14,24,25,26,27,28,29,30,31,32). We saw a very good correspondence between previous proposed diets for Mesozoic mammals and their position on the morphospace. See Supplementary Fig. 6 for a principal components analysis scatter plot which includes multituberculates and haramiyidans; these taxa were excluded from our study because the vast majority of them have jaw shapes dissimilar to the other extinct and extant mammals in our sample (i.e., allotherians have shorter jaws and thus more negative PC1 scores).
Stem mammals
Most stem mammals plot within the morphospace of extant insectivores and have positive PC1 scores. One exception is Sinoconodon (taxon #2, Fig. 3), which plots within the morphospace of extant carnivores; Sinoconodon is considered a carnivore based on dental morphology5. Haramiyavia (#1) is thought to have been a plant-dominated omnivore23 based on dental morphology, but here it plots within the morphospace of extant insectivores. Both morganucodontans in this study, Morganucodon (#3) and Dinnetherium (#4), have similar PC1 scores to extant insectivores, echoing the findings of Gill et al.14.
Molar morphology indicates omnivorous or faunivorous diets for docodontans; here they mostly plot within the morphospace of extant insectivores, with the exception of Haldanodon (#6) and Docofossor (#7). Agilodocodon (#9) was previously considered a plant-dominated omnivore, with exudativorous dental features which indicated a diet mainly composed of plant sap33; more recently, Wible and Burrows34 challenged this hypothesis and suggested that the teeth of Agilodocodon most closely resemble those of extant insectivores. Here, Agilodocodon plots firmly within the morphospace of extant insectivores, close to the insectivorous dusky antechinus (Antechinus swainsonii, #61) and the elephant shrews (Elephantulus rufescens [#114] and E. brachyrhynchus [#115]), which are insect-dominated omnivores.
According to Ji et al.28 the swimming docodontan, Castorocauda (#5), has dental features indicative of feeding on aquatic invertebrates and small vertebrates, like fish. Castorocauda is often depicted as being carnivorous and, particularly, piscivorous7,28,33. The jaw shape of Castorocauda is similar to that of modern day insectivores; this docodontan might have been feeding on “soft” aquatic invertebrates (Fig. 3). The other Mesozoic mammal in our sample proposed to have been semi-aquatic, Teinolophos (#13), plots in a similar area of the morphospace to Castorocauda. Our extant sample also includes a semi-aquatic carnivore, the water opossum (Chironectes minimus, #69), which plots in the middle of the carnivore morphospace, far away from Castorocauda and Teinolophos.
Docofossor (#7) has skeletal features indicative of a fossorial lifestyle and a dentition similar to those of extant mammals foraging underground, such as moles, solenodons, and tenrecs35. This docodontan has previously been considered an insectivore7. Here, Docofossor plots within the morphospace of extant carnivores; however, it plots close to the burrowing Hispaniolan solenodon (Solenodon paradoxus, #109), which has an insectivorous diet. Among the extant insectivores in our sample, the burrowing vermivores (e.g., the hairy-tailed mole, Parascalops breweri [#108], and the Hispaniolan solenodon) have more negative PC1 scores than other insectivores (similar to that of Docofossor), and their PC1 values are more similar to those of carnivores.
The dental morphology of Haldanodon (#6) is indicative of an insectivorous diet. Here, it plots within the carnivore morphospace (very near extant herbivores), because of its tall coronoid process and comparatively shorter jaw. Docodon (#8) likely ate insects and other small invertebrates27 and, based on its diminutive size36, Microdocodon (#10) was probably insectivorous. Both of these docodontans plot within the insectivore morphospace.
Non-therian crown mammals
The jaw shape of non-therian crown mammals varies widely, plotting mostly within the morphospace of insectivores and carnivores. Fruitafossor (#11), a fossorial mammal with teeth similar to extant armadillos, has been considered an omnivore eating insects, small invertebrates and some plants26. Here, it plots within the insectivore morphospace, closely to the insectivorous and fossorial hairy-tailed mole (Parascalops breweri, #108), and shares similar PC1 scores with other fossorial taxa, such as Docofossor (#7) and the Hispaniolan solenodon (#109).
Extant monotremes eat insects and other small invertebrates. It has been proposed that the Early Cretaceous monotreme Teinolophos (#13) had a semiaquatic lifestyle (on the basis of its enlarged mandibular canal37) and ate in a similar manner to the insectivorous Kuehneotherium38. Here Teinolophos, and the australosphenidan Henosferus (#12), have PC1 scores similar to insectivores and omnivores.
The eutriconodontans are a very diverse group of insectivores and carnivores which had a wide range of body sizes, including some of the largest Mesozoic mammals known1. Here all eutriconodontans fall within or very close to the extant carnivore morphospace. In particular, Triconodon (#16) and Argentoconodon (#19) plot within the carnivore morphospace, Trioracodon (#17) and Volaticotherium (#18) plot between the carnivore and insectivore morphospaces, and Yanoconodon (#15) plots within the insectivore morphospace. Both gobiconodontids, Gobiconodon (#20) and Repenomamus (#21), have more negative PC1 scores and plot closer to the herbivore morphospace, but still remain within or close to the carnivore morphospace. Triconodon, Trioracodon, Gobiconodon, and Repenomamus are all considered carnivores based on craniodental morphology and body size1,7,31; additionally, there is direct evidence for the carnivorous diet of Repenomamus from fossilized stomach contents4. Yanoconodon and Volaticotherium are considered insectivores7.
“Symmetrodontans” like Spalacotherium (#22), Zhangheotherium (#24) and Maotherium (#25) have often been considered insectivores based on their craniodental morphology1,7 (note “symmetrodontans” likely do not represent a monophyletic group, but are often grouped together based on their tooth morphology1). Here, all “symmetrodontans” plot within the insectivore morphospace. Dryolestids are also commonly considered insectivorous1,29. Here, Crusafontia (#26) plots between the morphospace of extant carnivores and insectivores, while Amblotherium (#27) plots within the insectivore morphospace.
Vincelestes (#29) has previously been considered a carnivore on the basis of jaw shape12. Here, it plots near the morphospaces of both omnivores and herbivores. Bonaparte24 considered the incisor wear of Vincelestes reminiscent of Cenozoic carnivores, and Rougier25 considered its jaw morphology indicative of a forceful bite enabling the incorporation of tough plant matter into a primarily carnivorous/insectivorous diet.
Therian crown mammals
Extant marsupials have a large diversity of diets, including herbivory, but the extinct metatherians in our sample are considered to have been limited in diet to insectivory and carnivory (note that there are some putatively herbivorous/omnivorous extinct metatherians, like Glasbius and polydolopimorphians39,40). Their jaw shape is very similar to that of extant carnivores and insectivores (Fig. 3). Dental morphology indicates that Eodelphis (#32) and Deltatheridium (#30) were carnivorous, Didelphodon (#31) durophagous or molluscivorous31,32, and Alphadon (#33) is considered to have been insectivorous, on the basis of its jaw shape and body size12. Dental microwear indicates a broad diet consisting of vertebrates, plants, and hard-shelled invertebrates for Didelphodon; biomechanical analyses of its skull and jaw points towards a durophagous diet15,16. Biomechanical analyses of the resistance to bending and torsion of Eodelphis jaws, points to a durophagous diet in Eodelphis cutleri and non-durophagous faunivory for Eodelphis browni16. Here, Eodelphis, Deltatheridium and Didelphodon plot closely to the extant carnivores, while Alphadon plots closely to the extant insectivores.
Extant placentals also have a wide range of diets, but many of the extinct eutherians in this study (i.e., Sinodelphys [#34], Juramaia [#35], Eomaia [#36], Kennalestes [#40], Barunlestes [#44], and Kulbeckia [#43]) are considered insectivorous7,12. Here, we corroborate this hypothesis (Fig. 3): all extinct eutherians plot within the insectivore morphospace, with the exception of Asioryctes (#38) which plots in the insectivore/carnivore morphospace, and Juramaia and Sinodelphys, which plot just outside the insectivore morphospace.
Using jaw shape to infer diet in Mesozoic mammals
We performed a phylogenetic flexible discriminant analysis (phylo FDA) following Motani and Schmitz41 to determine the posterior probability of the Mesozoic taxa belonging to one of three dietary categories: insectivore, carnivore, or herbivore (we omitted omnivores as they are not well discriminated in Fig. 3). We used the first seven PC scores (of the PCA of Procrustes coordinates of jaw shape), which together accounted for 81.39% of the variance. The results of the analysis can be seen in Fig. 4 and the posterior probability values can be seen in Supplementary Data 1. We used the extant taxa of known diets as the training dataset for the discriminant analysis: these taxa were classified correctly 89.19% of the time. For the most part, we see a good separation between dietary groups among extant mammals (Fig. 4a), with some exceptions: the primarily herbivorous olingo (Bassaricyon gabbii, #94) plots with the carnivores (although mainly frugivorous, it can consume small vertebrates), and a couple of insectivores plot very near the carnivores (i.e., the little brown bat [Myotis lucifugus, #104] and the Hispaniolan solenodon [Solenodon paradoxus, #109]). These three taxa, alongside the carnivorous greater bulldog bat (Noctilio leporinus, #101), were the only extant taxa misclassified by the discriminant analysis.
Fig. 4: Phylogenetic flexible discriminant analysis results, showing discriminant axis 1 (DA1) and two (DA2), of all taxa in this study.

Extinct taxa are color coded based on their posterior probability of belonging to one of the established dietary categories. Convex hulls show the position of the extant taxa in the plot and are color coded based on their dietary categories.

Full size image

The Mesozoic mammals included in our sample have largely been considered faunivorous and the results of the phylo FDA (Fig. 4b) corroborate this hypothesis. The majority of them are classified as insectivorous, including most stem mammals, australophenidans, “symmetrodontans” and eutherians, among others. Among the eutriconodontans, Argentoconodon, Gobiconodon, Repenomamus, and Trioracodon, are classified as carnivores, Triconodon and Yanoconodon are classified as insectivores, but with moderate support (posterior probabilities: 48% and 52%, respectively), and Phascolotherium and Volaticotherium are more confidently classified as insectivores (posterior probabilities: 60% and 73%, respectively). Among the metatherians, Didelphodon and Eodelphis are classified as carnivores, while Alphadon and Deltatheridium are classified as insectivores with moderate support (posterior probabilities: 54% and 52%, respectively). The stem mammals, Haramiyavia, Sinoconodon, and Docofossor are all confidently classified as carnivores (posterior probabilities over 80%), and the crown mammals Crusafontia and Kennalestes are also classified as carnivores, but with moderate support (posterior probabilities: 54% and 52%, respectively). Two taxa in the analysis are classified as herbivores, because of their relatively tall ascending rami: Vincelestes (#29) and Haldanodon (#6). The dental morphology of Vincelestes points to a primarily faunivorous diet24, but it has been previously noted that its jaw morphology is indicative of a forceful bite; Rougier25 suggested that this jaw morphology might have enabled Vincelestes to incorporate tough plant matter into its diet, but it might also be indicative of durophagy. The dental morphology27 and body size of Haldanodon point towards an insectivorous diet; in this analysis, the posterior probability of Haldanodon being a herbivore is not high (only 40.3%). The evidence thus far suggests Haldanodon had a faunivorous diet; its jaw morphology might be indicative of the incorporation of tougher food sources into its diet.
Mechanical advantage of the jaws of small mammals
We obtained mechanical advantage (MA) data to test whether extant mammals of different dietary groups have distinct MA values (Table 3). The mechanical advantage measurements were standardized across all jaws to account for differences in jaw morphology (e.g., presence or absence of the angular process) (Fig. 2b); the outlever was measured at the anterior end of the jaw and at the first lower molar (m1). When measuring mechanical advantage at the jaw tip and considering extant taxa only, we find statistically significant differences in the mechanical advantage of the masseter (MAM) values in all pairwise dietary combinations except for carnivore-insectivore (Table 3). The mechanical advantage of the temporalis (MAT) is statistically distinct only between herbivores and insectivores, and carnivores and insectivores (Table 3). Herbivores and carnivores do not have statistically distinct MAT values. This may differ in a sample of larger ( > 5 kg) therians. When measuring the outlever at the m1, we find statistically significant differences in all pairwise comparisons of MAM between dietary groups, except for herbivore–omnivore and carnivore–insectivore. When considering MAT, we only find significant differences between omnivores and carnivores, insectivores and herbivores, and insectivores and carnivores.
Table 3 Pairwise p values (uncorrected significance) of one-way PERMANOVAs of the mechanical advantage values of the masseter (MAM) and temporalis (MAT) obtained in this study on extant taxa of known dietary preferences only (permutation N = 9999).
Full size table

Figure 5 shows the mechanical advantage of the masseter (left) and temporalis (right), measured at the jaw tip, in a phylogenetic context (see also Supplementary Fig. 7 for individual taxon names). Phylogeny appears to have a large influence on the mechanical advantage and diet of the jaws of small mammals. Most Mesozoic taxa have low (blue) to intermediate (green) MAM values. Most stem mammals have intermediate (green) to high (red) MAM values and non-therian crown mammals have low MAM values, with the exception of Fruitafossor and Vincelestes (the latter has the highest MAM value of all taxa, both extinct and extant). Most eutherians, both extinct and extant, have intermediate to low MAM values, with the exception of the relatively high values (yellow to orange) seen in elephant shrews (order Macroscelidea) and the four-toed hedgehog (order Eulipotyphla, Atelerix albiventris). Some members of the orders Carnivora (including canids and euplerids) and Afrosoricida have some of the lowest MAM values. Metatherians have MAM values ranging from low to intermediate (in the orders Dasyuromorphia and Didelphimorphia, as well as in the Mesozoic metatherians) to some of the highest in the order Diprotodontia (e.g., the sugar glider [Petaurus breviceps], the woylie [Bettongia penicillata], the cuscus [Phalanger orientalis]).
Fig. 5: Mechanical advantage values of the masseter and temporalis when biting at the jaw tip visualized in the context of the phylogeny used in this study.

See Supplementary Fig. 7 for individual taxon names.

Full size image

Most taxa have intermediate MAT values (Fig. 5 and Supplementary Fig. 7). Very low MAT values are seen in the extinct non-therian crown mammals Teinolophos and Zhangheotherium and a few extant taxa, including marsupials like the Western barred bandicoot (Perameles bougainville) and the numbat (Myrmecobius fasciatus), and placentals such as the striped treeshrew (Tupaia dorsalis) and the short-snouted elephant shrew (Elephantulus brachyrhynchus). The highest MAT values belong to members of the order Carnivora, including skunks (Mephitis macroura and Conepatus humboldtii), the least weasel (Mustela nivalis) and the tayra (Eira barbara). Some diprotodontians like the common ringtail possum (Pseudocheirus peregrinus) and the sugar glider (Petaurus breviceps) also have relatively high MAT values. Some extinct taxa also have relatively high MAT values, including the stem mammal Docofossor, and the non-therian crown mammals, Triconodon and Vincelestes.
Figures 6 and 7 present a visualisation of the mechanical advantage of the masseter and the temporalis (x axis, outlever measured at the jaw tip) and the PC1 scores of Fig. 3 (y axis, mainly describes jaw length) because, as previously mentioned, this is the axis in which dietary categories among extant mammals are best discriminated. In the y axis of Figs. 6a and 7a, herbivores have short jaws, carnivores have short to intermediate-length jaws and insectivores have intermediate-length to long jaws. In Fig. 6a, insectivores and carnivores have low mechanical advantage values of the masseter (i.e., when biting: less forcefulness, more speed), and herbivores have higher mechanical advantage values (i.e., when biting: more forcefulness, less speed). In Fig. 7a, insectivores have lower mechanical advantage values of the temporalis, while carnivores and herbivores have higher mechanical advantage values. Note that most carnivores have intermediate MAT values, but some mustelids (i.e., the least weasel [Mustela nivalis, #99], the American badger [Taxidea taxus, #96], and the North American river otter [Lontra canadensis, #98]), have the highest MAT values among extant mammals. Also note that, among insectivores, those with the highest MAT values are burrowing vermivores (i.e., the short-tailed shrew tenrec [Microgale brevicaudata, #111], the hairy-tailed mole [Parascalops breweri, #108], and the Hispaniolan solenodon [Solenodon paradoxus, #109]). By using a combination of their MAM and MAT values (as well as their jaw length), we can distinguish dietary categories among extant mammals. We decided to omit omnivores from these figures because, as seen in Fig. 3, they cannot be distinguished from other dietary groups on the basis of jaw shape.
Fig. 6: Scatter plot of the mechanical advantage of the masseter (x axis) vs. PC1 scores from Fig. 3 (y axis), which mainly describes jaw length.

a Extant taxa, b extinct taxa. Colors indicate known dietary categories of extant mammals and suggested dietary categories for Mesozoic mammals (obtained from the literature). Ovals indicate where extant taxa of known dietary categories plot, as in part a.

Full size image

Fig. 7: Scatter plot of the mechanical advantage of the temporalis (x axis) vs. PC1 scores from Fig. 3 (y axis), which mainly describes jaw length.

a Extant taxa, b extinct taxa. Colors indicate known dietary categories of extant mammals and suggested dietary categories for Mesozoic mammals (obtained from the literature). Ovals indicate where extant taxa of known dietary categories plot, as in part a.

Full size image

We also obtained additional mechanical advantage measurements, in which the outlever was measured at the first lower molar (m1), rather than the jaw tip (Supplementary Figs. 8, 10, 11, and 13). We made this alternative measurement because Grossnickle17 found that the length of the posterior portion of the jaw (measured from the jaw joint to the m1) is a strong predictor of diet in mammals. Compared to the mechanical advantage (MA) measurements at the jaw tip (Figs. 6a and 7a), we see a less clear distinction between dietary groups among extant mammals. There is considerable overlap between dietary groups in Supplementary Fig. 10 (jaw length~MAM). In Supplementary Fig. 11 (jaw length~MAT), there is a better separation between dietary groups.
Based on previous likely determinations of diet of Mesozoic mammals (see Supplementary Data 1 for the full list of sources), most taxa plot where it is “expected” of them, with some exceptions (Figs. 6b and 7b): 1) about half of the stem mammals (i.e., Haramiyavia, Sinoconodon, Morganucodon, Haldanodon, and Docofossor), most of which are thought to have been faunivorous, have higher MAM values than modern insectivores and carnivores, and 2) the docodontan Castorocauda has MAM and MAT values consistent with an insectivorous diet, as opposed to the carnivorous diet proposed for this taxon7,28,33. Most Mesozoic mammals have mechanical advantage values similar to modern insectivores, a few taxa are similar to carnivores (e.g., Sinoconodon, Triconodon, Trioracodon, Argentoconodon, Gobiconodon, Repenomamus, Deltatheridium, Didelphodon, and Eodelphis), and some are more similar to herbivores (e.g., Vincelestes and Fruitafossor). More

Model equations
We introduce a spatiotemporal model to describe the bacteria-phage interaction in the upper part of the soil with the depth H (we consider (H=1) m) in a typical agricultural field. Here we consider a 1D model where all abiotic and biotic components depend on time t and vertical coordinate h. The biotic component of the model consists of 4 compartments: phage-free bacteria (S) susceptible to infection by the phage, bacteria infected by the phage in its lysogenic ((I_1)) and lytic ((I_2)) states, and free phages (P). The total density of the host bacterial populations N is defined as (N = S + I_1 + I_2). The schematic diagram illustrating bacteria-phage interactions is similar to that of Egilmez and co-authors16. The local species interactions are described based on the classical modelling approach6,19. Our spatiotemporal model is of reaction-diffusion type and is described by the following equations

$$begin{aligned} begin{aligned} frac{partial S(t,h)}{partial t}&= D_b frac{partial ^2 S(t,h)}{partial h^2} +alpha (T) S(t,h) Big [1-frac{N(t,h)}{C(h)}Big ] – K_S S(t,h)P(t,h), \ frac{partial I_1(t,h)}{partial t}&= D_b frac{partial ^2 I_1(t,h)}{partial h^2} + {overline{alpha }}(T) I_1(t,h) Big [1- frac{N(t,h)}{C(h)}Big ] + K_1(T) S(t,h) P(t,h) – lambda _1(T) I_1(t,h), \ frac{partial I_2(t,h)}{partial t}&= D_b frac{partial ^2 I_2(t,h)}{partial h^2} + K_2(T) S(t,h) P(t,h) + lambda _1(T) I_1(t,h) – lambda _2 I_2(t,h), \ frac{partial P(t,h)}{partial t}&= D_P frac{partial ^2 P(t,h)}{partial h^2} -K N(t,h) P(t,h) – mu P(t,h) + b lambda _2 I_2(t,h). end{aligned} end{aligned}$$
(1)

In the above model, we parameterise the growth of susceptible bacteria via a standard logistic growth function6, where (alpha) is the maximal per capita growth rate and C is the carrying capacity of the environment; we assume that C(h) varies with depth. Infection of S by phages P at low temperatures results in lysogeny which is described by a mass action term (K_s S(t,h) P(t,h)). The growth of lysogenic bacteria (I_1) is described by a logistic function as in the case of S; however, with a different maximal growth rate ({overline{alpha }} (T)) as detailed in the next subsection. At high temperatures, the transition from the lysogenic to the lytic cycle of infection occurs: this is described by the term (lambda _1 (T) I_1(t,h)). Infection by the phage via the lytic cycle is modelled by the term (K_2 (T)S(t)P(t)). The death rate of infected bacteria due to lysis is modelled by (lambda _2 (T) I_2). The lysis of a bacterium results in the release of b new phages, the the burst size6. In the equation for P, KN(t)P(t) stands for the loss of phage due to binding to any type of bacteria (for simplicity, we assume that there is no saturation in binding). Finally, (mu P(t,h)) is the natural mortality or deactivation of phages.
According to this framework, the vertical displacement of the phage and bacteria are modelled by a diffusion term (first term in each equation), where (D_b) and (D_P) are the diffusion coefficients of bacteria and phage, respectively. The variation of the temperature T across the soil is described by the heat equation

$$begin{aligned} frac{partial T(t,h)}{partial t} = D_h frac{partial ^2 T(t,h)}{partial h^2}, end{aligned}$$
(2)

where (D_h) is the diffusion coefficient of heat transfer (see more detail in the next section). Models (1)–(2) should be supplied with appropriate boundary conditions. We assume that the model has the zero-flux boundary condition for all biotic components (bacteria and phage) at (h=0) and (h=H). For the temperature, we consider Dirichlet boundary conditions such that (T(t,0)= T_s (t)) and (T(t,H)= T_H), where (T_s (t)) is the surface temperature and (T_H) is a constant temperature in deeper soil layers.
Parameterisation of equation terms
Next we describe the functional forms of the dependence of model parameters on the temperature and the depth. Following the previous study16, we assume that the maximal bacterial growth rates (alpha (T)) and ({overline{alpha }}(T)) are described by

$$begin{aligned} alpha (T)= & {} exp left (-frac{(T-T_0)^2}{2sigma ^2}right )alpha _{text {max}}, end{aligned}$$
(3)

$$begin{aligned} {overline{alpha }}(T)= & {} alpha (T) left [1-frac{T^n}{T_1^n + T^n}right ] = alpha _{text {max}}exp {left (-frac{(T-T_0)^2}{2sigma ^2}right )} left [1-frac{T^n}{T_1^n + T^n}right ], end{aligned}$$
(4)

where (T_0=38.2 ^circ text {C}) is the optimal temperature; (T_1=34.8 ^circ text {C}) is the temperature corresponding to the switch between the lytic and the lysogenic cycles; (alpha _{text{max}}=23 text {day}^{-1}) is the maximal possible growth, (sigma =9.1 ^circ text {C}) describes the decay of growth with temperature T16,20.
In the equation for ({overline{alpha }}(T)), we assume that at a high temperature normal cell division of (I_1) stops since there is a transition to a lytic state in bacteria. In the soil bacteria grow anaerobically or microaerophillically, and the growth rates of B. pseudomallei under such conditions are yet to be studied. For simplicity they are assumed to be the same as under aerobic conditions. Realistic values of the above parameters are listed in Table 1. Note that in the model both (alpha (T)) and ({overline{alpha }}(T)) are in fact effective growth rates of the bacterial populations, i.e. they incorporate the replication of cells and as well as their mortality.
Table 1 Parameters used in the model along with their units and ranges.
Full size table

The overall adsorption rate of the phage K is estimated as (2 times 10^{-7} text {ml}^{-1} text { day}^{-1}) from Egilmez et al.16. The adsorption constants (K_1 (T)), (K_2 (T)) and the transition rate from lysogenic to lytic cycle (lambda _1(T)) depend on temperature as follows16:

$$begin{aligned} K_1(T)&= left(1-frac{T^n}{T_1^n + T^n}right) K_S, nonumber \ K_2(T)&= frac{T^n}{T_1^n + T^n} K_S, nonumber \ lambda _1(T)&= frac{T^n}{T_1^n + T^n} {lambda _1}_{text {max}} , end{aligned}$$
(5)

where (K_S) is the maximal phage adsorption constant ((K_S=epsilon K) where (epsilon =0.3) is the adsorption efficiency) and (lambda _{1_text {max}}=23 text {day}^{-1}) is the maximal transition rate which is assumed to be equal to the maximal growth rate of the bacteria16. The lysis rate of bacteria (lambda _2=20 text {day}^{-1}) (depending on 50 min latency time13) and the burst size (b = 100) in the model are assumed to be constant16. The temperature dependence of (alpha (T)), ({overline{alpha }}(T)), (K_1 (T)), (K_2 (T)) and (lambda _1(T)) are shown in Fig. 2. The mortality rate of phages (mu) is high near the surface due to ultraviolet radiation, but the role of ultraviolet radiation becomes negligible starting from a depth of a few centimetres because sunlight cannot penetrate the soil. For the above reason, we can assume (mu =3 text {day}^{-1}) to be constant.
Figure 2

(a) Temperature dependence of the adsorption constants (K_i) ((i=1,2)) of the phage (measured in (text {ml}^{-1} text {day}^{-1})). (b) Growth rates of susceptible (alpha (T)) and lysogenic ({overline{alpha }}(T)) bacteria and the transition rate (lambda _1(T)) from the lysogenic cycle to the lytic cycle (measured in (text {day}^{-1})). The corresponding analytical expressions for the temperature dependence are given by (3)–(5).

Full size image

The carrying capacity C of the bacteria varies with the depth of the soil, according to empirical observations21,22,23. This can be explained by the fact that the humus, oxygen, nitrogen contents, or/and water content in the soil generally decrease with depth24. We use a combined approach to parameterise C(h) based on the available empirical data. We assume that in the absence of phages, the bacteria achieve numbers close to the carrying capacity at a given depth. Firstly, we parameterise the dependence of the overall bacterial load on depth in paddy soils in Southern Asia using the existing data22. Then we re-scale the obtained curve based on the available observations of B. pseudomallei at a depth (h=30 text {cm})25,26. We approximate C(h) using the following simple Gaussian-type curve

$$begin{aligned} C(h)=(C_text {surf} -C_0)exp (-B h^2)+C_0, end{aligned}$$
(6)

where (C_text {surf}) gives the maximal number of bacteria near the surface (h), B determines how fast the bacterial abundance decreases with depth, (C_0) is background bacterial density which takes into account the fact that bacteria can survive even at large depths (e.g. (h=100 text {cm})). Based on our estimates (see supplementary material SM1 for more detail), we will use the following parameter values as defaults: (C_text {surf} = 1 times 10^6) (text {cell/ml}), (B=7.5 times 10^{-4}) (1/{text{cm}}^2), (C_0=10^4) (text {cell/ml}) . One can easily see that C(h) has a maximum at the surface and monotonically decreases with depth. We assume that the carrying capacity of the environment is not influenced by seasonal variation.
The coefficient (D_h) in the equation for the temperature distribution can be estimated as follows. Generally, (D_h) is related to (rho _s), (C_{rho s}) and (k_s) which are the bulk density, specific heat and thermal conductivity in soil, respectively, i.e. (D_h=k_s/(rho _s C_{rho s)}). We use the estimates for (rho _s), (C_{rho s}) and (k_s) from27 which gives (rho _s=110.52 text {kg}/text{m}^3), (C_{rho s} = 1130) (text {J/kg K}) and (k_s = 0.0967) (text {W/m K}) and, for the diffusion coefficient (D_h=7.7 times 10^{-8}) (text {m}^2 text{s}^{-1}). The variation of (T_s)—the surface temperature—is obtained from the historical weather report for the surface16. The bottom boundary temperature (T_H) at (h=H=1 text {m}) is considered to be (22 ^circ text{C}). The initial value of the temperature distribution (T_s (0)) is assumed to be linear, but this assumption does not affect long-term temperature dynamics.
The paddy fields in which we model the bacteria-phage interactions are flooded lands, where the soil is either mud or muddy water. Many factors can affect vertical dispersal of bacteria and phages in such soil. For instance, rain water can carry bacteria and phage up or down in the soil, which can be mathematically modelled by adding an advection term; however, for simplicity we ignore such effects in this paper. We also assume the phage and bacteria vertical diffusion coefficients to be constant; however, it is rather hard to provide accurate estimates for (D_p) and (D_b). In water, the diffusion coefficient of bacteria and phages can be estimated as (3.6times 10 ^{-10} text {m}^2 text{s}^{-1}= 0.3 text {cm}^2 text{day}^{-1}) and (2.8 times 10^{-12} text {m}^2 text{s}^{-1}= 0.002 text {cm}^2 text{day}^{-1}), respectively28, but the diffusivity in soil should be smaller than this. As such, these values should be considered as upper limits for (D_P) and (D_b), with the actual coefficients being orders of magnitude smaller. We undertook simulations with different combinations of diffusion coefficients in this range, and found that the patterns of vertical distribution do not largely depend on the diffusion coefficients provided (D_P< 10^{-3} text {cm}^2 text{day}^{-1}) and (D_b < 10^{-2} text {cm}^2 text{day}^{-1}), due to the strong external forcing of the system by temperature (see “Results” section for details). In our numerical simulations, we use both explicit and implicit numerical schemes. We take a 0.1 cm spatial step size to get a proper resolution. We separately compute the heat equation to define T(t) with a smaller time resolution and then apply the temperature obtained to model bacteria-phage interactions for a larger time resolution (for example, (Delta t cong 7 times 10^{-5}) day). We compute the average densities of the species (both in terms of spatial and temporal averaging) using a numerical right Riemann sum. The accuracy of our numerical simulation was verified by reducing both time and space steps and comparing the results obtained. We use daily and seasonal variation of temperatures (for the period of 2013–2016) in the provinces of Nakhon Phanom and Sa Kaeo in Thailand to parameterise the model (http://www.worldweatheronline.com). The unit of the densities of bacteria and phages are cells/ml. The summary of model parameters as well their values are provided in Table 1. More

Elements of TA
Identifying biogeographic trends within volatile data required several innovations in the MASTIF model20, building from multivariate state-space methods in previous applications41,52. Standard modeling options, such as generalized linear models and their derivatives, do not accommodate key features of the masting processes. First, multiple data types are not independent. Maturation status is binary with detection error, CCs are non-negative integers, also with detection error, and STs require a transport model (dispersal) linking traps to trees, and identification error in seed identification. Of course, a tree observed to bear seed, now or in the past, is known to be mature now. However, failure to observe seed does not mean that an individual is immature because there are detection errors and failed crop years41,64.
Second, seed production is quasiperiodic within an individual (serial dependence), quasi-synchronous between individuals (“mast years”), [and] there is dependence on environmental variation, and massive variation within and between trees41,53,65. Autoregressive error structures (AR(p) for p lag terms) impose a rigid assumption of dependence that is not consistent with quasiperiodic variation that can drift between dominant cycles within the same individual over time43. It does not allow for individual differences in mast periodicity.
Third, climate variables that affect fecundity operate both through interannual anomalies over time and as [a] geographic variation. The masting literature deals almost exclusively with the former, but our application must identify the latter: the potentially smooth variation of climate effects across regions must be extracted from the many individual time series, each dominated by local “noise.”
Finally, model fitting is controlled by the size classes that dominate a given site and thus is insensitive to size classes that are poorly represented. Large trees are relatively rare in eastern forests, making it hard to identify potential declines in large, old individuals41,53. Conversely, the shade-intolerant species that dominate second-growth forests often lack the smaller size classes needed to estimate maturation and early stages where fecundity may be increasing rapidly.
Several of the foregoing challenges are resolved in the MASTIF model by introducing latent states for individual maturation status and tree-year seed production. The dependent data types (maturation status, CCs, STs) become conditionally independent in the hierarchical MASTIF model (e.g., ref. 66). The serial dependence is handled as a conditional hidden Markov process for maturation that combines with CCs and STs by way of stochastic (latent) conditional fecundity. Maturation status and conditional fecundity must be estimated jointly, that is, not with separate models. The latent maturation/fecundity treatment avoids imposing a specific AR(p) structure. In the MASTIF model there is a posterior covariance in maturation/fecundity across all tree-year estimates that need not adhere to any specific assumption20. The dependence across individuals and years is automatic and available from the posterior distribution.
Separating the spatial from temporal components of climate effects is possible here, not only because the entire network is analyzed together but also because predictors in the model include both climate norms for the individual sites and interannual anomalies across sites35,52. TA depends on both of these components.
Extracting the trends from volatile data further benefits from random individual effects for each tree and the combination of climate anomalies and year effects over time. A substantial literature focuses on specific combinations of climate variables that best explain year-to-year fecundity variation, including combinations of temperature, moisture, and water balance during specific seasons over current and previous years19,25,41. Results vary for each study, presumably due to the differences in sites, species, size classes, duration, data type, and modeling assumptions. For TA, the goal is to accommodate the local interannual variation to optimize identification of trends in space and time. Thus, we include a small selection of important climate anomalies (spring minimum T of the current year, summer T of the current and previous year, and moisture D of the current and previous year). The climate anomalies considered here do not include every variable combination that could be important for all size classes of every species on every site. For this reason, we combine climate anomalies with year effects. Year effects in the model are fixed effects within an ecoregion and random between ecoregions (ecoregions are shown in Fig. 2 and listed in Supplementary Data 2). They are fixed within an ecoregion because they are not interpreted as exchangeable and drawn at random from a large population of possible years. They are random between ecoregions due to the uneven distribution of sites (Supplementary Data 1)20.
To optimize inference on size effects, the sampling of coefficients in posterior simulation is implemented as a weighted regression. This means that the contribution of tree diameter to fecundity is inversely proportional to the abundance of that size class in the data. This approach has the effect of balancing the contributions of abundant and rare sizes. Identifying size effects further benefits from the introduction of opportunistic field sampling, which can target the large individuals that are typically absent from field study plots.
MASTIF data network
Data included in the analysis come from published and unpublished sources and offer one or both of the two data types, CCs and STs (Supplementary Data 1). Both data types inform tree-year fecundity; they are plotted by year in Fig. 6.
Fig. 6: Distribution of observation trees by year in the North American region of the MASTIF network.

Sites are listed by ecoregion in the Supplementary Data 2.

Full size image

CCs in the MASTIF network are obtained by one of three methods. Most common are counts with binoculars that are recorded with an estimate of the fraction of the crop that was observed. A second CC method makes use of seeds collected per ground surface area relative to the crown area. This method is used where conspecific crowns are isolated and wind dispersal is limited. The crop fraction is the ratio of ground area for traps relative to the projected crown area. Examples include HNHR67 and BCEF68.
A third CC method is based on evidence for past cone production that is preserved on trees. This has been used for Abies balsamea at western Quebec sites69, Pinus ponderosa in the Rocky Mountains70, and for Pinus edulis at SW sites27.
ST data include observations on individual trees that combine with seed counts from traps. Because individual studies can report different subcategories of seeds, and few conduct rigorous tests of viability, we had to combine them using the closest description to the concept of “viable”. For example, we do not include empty conifer seeds. A dispersion model provides estimates of seeds derived from each tree. ST and CC studies are listed in Supplementary Data 1. The likelihoods for CCs and STs are detailed in ref. 20. Individually and in combination, the two data types provide estimates, with full uncertainty, for fecundity across all tree-years.
Fitted species had multiple years of observations from multiple sites, which included 211,146 trees and 2,566,594 tree-years from 123 species. Sites are shown in Fig. 2 of the main text by ecoregion, they are named in Fig. 1 and summarized in Supplementary Data 1. For TA the fits were applied to 7,723,671 trees on inventory plots. Mean estimates for the genus were used for inventory trees belonging to species for which there were not confident fits in the MASTIF model, which amounted to 7.2% of inventory trees. Detailed site information is available at the website MASTIF.
Covariates
Covariates in the model include as main effects tree diameter, tree canopy class (shading), and the climate variables in Fig. 1 of the main text and described in Table 1. A quadratic diameter term in the MASTIF model allows for changes in diameter response with size52. Shade classes follow the USDA Forest Inventory and Analysis (FIA)/National Ecological Observation Network (NEON) scheme that ranges from a fully exposed canopy that does not interact with canopies of other trees to fully shaded in the understory. Shading provides information on competition that has proved highly significant for fecundity in previous analyses41,52.
Table 1 Predictors in the model, not all of which are important for all species.
Full size table

To distinguish between the effects of spatial variation versus interannual variability, spring T and moisture D are included in the model as site means and site anomalies35. Spring minimum T affect phenology and frost risk during flowering and early fruit initiation. Summer mean T (June–August) is included both as a linear and quadratic term. Mean summer T is linked to thermal energy availability during the growing season, with the quadratic term allowing for potential suppression due to extreme heat. Moisture D (cumulative monthly PET-P (potential evapotranspiration[-] minus precipitation) for January–August) is included as a site mean and an annual anomaly. Moisture D is important for carbon assimilation and fruit development during summer in the eastern continent and, additionally, from the preceding winter in the western continent. For species that develop over spring and summer, anomalies incorporate the current and previous year. We did not include longer lags in covariates. For species that disperse seed in spring (Ulmus spp. and some members of Acer), only the previous year was used. Temperature anomalies were included for spring, but not summer, simply to reduce the number of times that temperature variables enter the model, and these two variables tended to be correlated at many sites.
Climate covariates were derived from gridded climate products and combined with local climate monitoring where it is available. Terraclimate71 provides monthly resolution, but it is spatially coarse. For both norms and trends, we used the period from 1990 to 2019 because global temperatures have been increasing consistently since the 1980s, and this span broadly overlaps with fecundity data (Fig. 6). CHELSA72 data are downscaled to a 1 km grid, but it does not extend to 2019. Our three-component climate scaling used regression to project CHELSA forward using Terraclimate, followed by downscaling to 1 km with CHELSA, with further downscaling to local climate data. Even where local climate data exist, they often do not span the full duration of field studies, making the link to gridded climate data important. Local climate data were especially important for mountainous sites in the Appalachians, Rockies, Sierra Nevada, and Cascades.
Of the full list of variables, a subset was retained, depending on species (some have narrow geographic ranges) and deviance information criteria of the fitted model (Supplementary Data 2). Year effects in the model allow for the interannual variation that is not absorbed by anomalies20.
Model fitting and TA
As mentioned above, model fitting applied the hierarchical Bayes model of ref. 20 to the combination of time series and opportunistic observations summarized in Fig. 1. Posterior simulation was completed with Markov chain Monte Carlo based on direct sampling, Metropolis, and Hamiltonian Markov chain. Model fitting used 211,146 trees and 2,566,594 tree-years from 123 species (Supplementary Data 2). Only species with multiple observation years were included.
The climate variable referenced as C in Eq. (1) of the main text is, in fact, a vector of climate variables described in the previous section, spring minimum T, summer mean T, and moisture D (Table 1). The anomalies and year effects in the fitted model contribute to the trends not explained by biogeographic variation as γ in Eq. (1). For main effects in the model, the partial derivatives are fitted coefficients, an example being the response to spring minimum temperature (partial f/partial {T}_{mathrm{sp}}={beta }_{{T}_{mathrm{sp}}}). For predictors involved in interactions, the partial derivatives are combinations of fitted coefficients and variables. For example, the response to moisture D, which interacts with tree size, is (partial [F], f/partial {D}={beta }_{{D}} + beta_{GD}G). The response to diameter G, which is quadratic and interacts with D, is (partial f/partial G={beta }_{G}+2{beta }_{{G}^{2}}G ,+{beta }_{GD}D).
Trend decomposition applied the fitted model to every tree in forest inventories from the USDA FIA program, the Canada’s National Forest Inventory, the NEON, and our MASTIF collaboration. Each tree in these inventories has a species and diameter. For trees that lack a canopy class, regression was used to predict it from distances and tree diameters based on inventories that include both location and canopy class, including NEON, FIA, and the MASTIF network. Although inventories differ in the minimum diameter they record, few trees are reproductive at diameters below the lower diameter limits in these surveys, so the effect on fecundity estimates is negligible. For the indirect effects of climate coming through tree growth rates, the same covariates were fitted to growth as previously defined for fecundity, using the change in diameter observed over multiple inventories. A Tobit model was used to accommodate the fact that a second measurement can be smaller than an earlier measurement. The Tobit thus treats negative growth as censored at zero. TA to inventory plots used 7,717,677 trees. Because not all species in the inventory data are included in the MASTIF network, mean fecundity parameters for the genus were used for unfitted species. Species fitted in the MASTIF network accounted for >90% of trees in inventory plots (Supplementary Data 2).
From the predictive distributions for every tree in the inventory data, we evaluated predictive mean trends aggregated to species and plot in Fig. 2b. We extracted specific terms that comprise the components in Fig. 4 and aggregated them too to the plot averages.
General form for TA
Equation 1 simplifies the model to highlight direct and indirect effects. Again, climate variables and tree size represent only a subset of the predictors in the model that are collected in a design vector ({{bf{x}}}_{t}=[{x}_{1,t},ldots ,{x}_{Q,t}]^{prime}), where the q = 1, …, Q predictors include shading from local competition, individual size, and climate and habitat variables (Table 1). On the proportionate scale, Eq. (1) can be written in terms of all predictors, including main effects and interactions, as

$$frac{{mathrm{d}}f}{{mathrm{d}}t}=mathop{sum }limits_{q=1}^{Q}left(frac{partial f}{partial {x}_{q}}+sum _{q^{prime} in {I}_{q}}frac{partial f}{partial ({x}_{q}{x}_{q^{prime} })}{x}_{q^{prime} }right)frac{{mathrm{d}}{x}_{q}}{{mathrm{d}}t}+gamma$$
(2)

where Iq are variables that interact with xq. In this application, interactions include tree diameter with moisture deficit and diameter squared. Each term in the summation consists of a main effect of xq and interactions that are multiplied by the rate of change in variable xq. For the simple case of only two predictors, Eq. (2) is recognizable as Eq. (1) of the main text, where x1, x2 have been substituted for variables G and C. In our application, predictors include additional climate and shading (Table 1).
Recognizing that environmental variables affect not only fecundity but also growth rate, we extract the size effect, that is, the xq that is G, and incorporate these indirect effects (through growth) by expanding g = dG/dt in Eq. (1) of the main text as

$$g=mathop{sum }limits_{q=1}^{Q}left(frac{partial g}{partial {x}_{q}}+mathop{sum}limits _{q^{prime} in {I}_{q}}frac{partial g}{partial ({x}_{q}{x}_{q^{prime} })}{x}_{q^{prime} }right){x}_{q}+nu$$
(3)

where ν is the component of growth that is not accommodated by other terms. This expression allows us to evaluate the full effect of climate variables, including those coming indirectly through growth.
Connecting fitted coefficients in MASTIF to TA
This section connects the continuous, deterministic Eq. (1) to the MASTIF model20 with the interpretation of responses, direct effects, and full effects of Fig. 5. To summarize key elements of the fitted model20, consider a tree i at site j that grows to reproductive maturity and then produces seed depending on its size, local competitive environment, and climate. We wish to estimate the effects of its changing environment and condition on fecundity using a model that includes spatial variation in predictors that are tracked longitudinally over years t. Fecundity changes through maturation probability ρij(t), which increases as trees increase in size, and through conditional fecundity ψij(t), the annual seed production of a mature tree. Let zij(t) = 1 be the event that a randomly selected tree i is mature in year t. Then, ρij(t) is the corresponding probability that the tree is mature, E[zij(t)] = ρij(t)(ρ is not to be confused with the probability that a tree that is now immature will make the transition to the mature state in an interval dt = 1. That is a different quantity detailed in the Supplement to ref. 41). Fecundity has expected value Fij(t) = ρij(t)ψij(t). On a proportionate (log) scale,

$${f}_{ij}(t)={mathrm{log}},{F}_{ij}(t)={mathrm{log}},{rho }_{ij}(t)+{mathrm{log}},{psi }_{it}(t)$$
(4)

The corresponding rate equation is

$$frac{{mathrm{d}}f}{{mathrm{d}}t}=frac{{mathrm{d}},{mathrm{log}},rho }{{mathrm{d}}t}+frac{{mathrm{d}},{mathrm{log}},psi }{{mathrm{d}}t}$$
(5)

The discretized and stochasticized version of Eq. (1) is

$$frac{{mathrm{d}}{F}_{ij}}{{mathrm{d}}t} = , frac{{F}_{ij,t+{mathrm{d}}t}-{F}_{ij,t}}{{mathrm{d}}t}+{epsilon }_{ij,t}\ = , {{Delta }}{F}_{ij,t}+{epsilon }_{ij,t}$$
(6)

where dt = 1 and ϵij,t is the integration error. When applied to a dynamic process model, this term further absorbs process error (see above), which is critical here to allow for conditional independence where observations are serially dependent. In simplest terms, ϵ is model miss-specification that allows for dependence in data.
The MASTIF model that provides estimates for TA is detailed in ref. 20. Elements of central interest for TA are

$${F}_{ij,t} = , {z}_{ij,t}{psi }_{ij,t}\ left[{z}_{ij,t}=1right] sim , {{Bernoulli}}left({rho }_{ij,t}right)\ {rho }_{ij,t} = , {{Phi }}({{boldsymbol{mu }}}_{ij,t})\ mathrm{log},{psi }_{ij,t} = ,{{bf{x}}}_{ij,t}^{prime}{boldsymbol{beta }}+{h}_{t}left(Tright)+{epsilon }_{ij,t}$$

where μij,t = α0 + αGGij,t describes the increase in maturation probability with size, Φ(⋅) is the standard normal distribution function (a probit), ϵij,t ~ N(0, σ2), and ht(T) can include year effects, h(T) = κt, or lagged effects, (h(T)=mathop{sum }nolimits_{k = 1}^{p}{kappa }_{k}{psi }_{ij,t-k}), that contribute to γ in Eq. (1) of the main text. If year effects are used, then γ includes the trend in year effects. (The generative version of this model writes individual states at t conditional on t − 1 and is given in the Supplement to ref. 20.). If an AR(p) model is used, then γ = κ1 (provided data are not detrended). Random individual effects in the fitted model are marginalized for prediction of trees that were not fitted, meaning that σ2 is the sum of model residual and random-effects variance. Again, the length-Q design vector xij,t includes individual attributes (e.g., diameter Gij,t), local competitive environment, and climate (Table 1). There is a corresponding coefficient vector β.
Moving to a difference equation (rate of change) for conditional log fecundity,

$${{Delta }}{f}_{ij,t}={{Delta }}mathrm{log},{rho }_{ij,t}+{{Delta }}mathrm{log},{psi }_{ij,t}$$

where

$${{Delta }}mathrm{log},{psi }_{ij,t} =mathrm{log},{psi }_{ij,t+1}-mathrm{log},{psi }_{ij,t}\ ={{Delta }}{{bf{x}}}_{ij,t}^{prime}{boldsymbol{beta }}+{gamma }_{ij,t}+{nu }_{ij,t}\ {{Delta }}{{bf{x}}}_{ij,t} ={{bf{x}}}_{i,t}-{{bf{x}}}_{ij,t-1}\ {nu }_{ij,t} sim N(0,2{sigma }^{2})$$

The variance in the last line is the variance of the difference Δϵij,t.
Elements of basic theory in Eq. (1) of the main text are linked to data through the modeling framework as

$${{Delta }}{f}_{ij,t}= +{beta }_{{T}_{sp}}{{Delta }}{T}_{sp,j}\ +left({beta }_{T}+2{beta }_{{T}^{2}}{T}_{j}right){{Delta }}{T}_{j}\ +left({beta }_{D}+{beta }_{GD}{G}_{ij,t}right){{Delta }}{D}_{j}\ +left({alpha }_{G}frac{phi ({{boldsymbol{mu }}}_{ij,t})}{{{Phi }}({{boldsymbol{mu }}}_{ij,t})}+{beta }_{G}+2{beta }_{{G}^{2}}{G}_{ij,t}+{beta }_{GD}{D}_{j}right){{Delta }}{G}_{ij}\ +{gamma }_{ij,t}+{nu }_{ij,t}$$
(7)

where ϕ(⋅) is the standard normal density function that comes from the rate of progress toward maturation. Again, the anomalies do not appear in this expression for trends because trends in the anomalies and year effects enter through γ.
The first four lines in Eq. (7) are, respectively, the effects of trends in spring minimum temperatures ΔTsp,j, summer mean temperature ΔTj, moisture deficit ΔDj, and size ΔGij, where the latter comes from growth on inventory plots. The contribution of maturation to change in fecundity is the first term in the fourth line, αGϕ(μij,t)/Φ(μij,t). A map of this term in Fig. 7b shows the strong contribution to fecundity in the East due to the young (Fig. 7a) and/or small (Fig. 4b) trees there. The sum of these terms dominates the patterns in Fig. 3c.
Fig. 7: Size and maturation effects on fecundity.

a Stand age variable in FIA data and b positive effect of maturation for increasing fecundity in the eastern continent. In the West, maturation does not contribute to rising fecundity because large trees are predominantly [mature] larger.

Full size image

All terms in Eq. (7) have units of mean change in proportionate fecundity, and these are mapped in figures of the main text. We focus on proportionate fecundity because it reflects the full effect of climate as opposed to total fecundity, which would often be dominated by one or a few trees of a single species. However, from proportionate fecundity we can obtain change in fecundity as ΔFij,t = Fij,t × Δfij. Stand-level effects on fecundity change at site j can be obtained from individual change as

$${{Delta }}{F}_{j}=mathop{sum }limits_{i=1}^{{n}_{j}}{{Delta }}{F}_{ij}=mathop{sum }limits_{i=1}^{{n}_{j}}{F}_{ij}{{Delta }}{f}_{ij,t}$$

Again, maps in Fig. 5 show mean proportionate effects for all trees on an inventory plot.
Reporting summary
Further information on research design is available in the Nature Research Reporting Summary linked to this article. More

1.
Ockendon, N. et al. Mechanisms underpinning climatic impacts on natural populations: altered species interactions are more important than direct effects. Glob. Chang Biol. 20, 2221–2229 (2014).
ADS  PubMed  Article  PubMed Central  Google Scholar 
2.
Boggs, C. L. & Inouye, D. W. A single climate driver has direct and indirect effects on insect population dynamics. Ecol. Lett. 15, 502–508 (2012).
PubMed  Article  PubMed Central  Google Scholar 

3.
Burkett, V. R. et al. Nonlinear dynamics in ecosystem response to climatic change: case studies and policy implications. Ecol. Complex. 2, 357–394 (2005).
Article  Google Scholar 

4.
Molnár, P. K., Sckrabulis, J. P., Altman, K. A. & Raffel, T. R. Thermal performance curves and the metabolic theory of ecology—a practical guide to models and experiments for parasitologists. J. Parasitol. 103, 423–439 (2017).

5.
Hortion, J. et al. Acute flavivirus and alphavirus infections among children in two different areas of Kenya, 2015. Am. J. Trop. Med. Hyg. 100, 170–173 (2019).
PubMed  Article  PubMed Central  Google Scholar 

6.
Stewart-Ibarra, A. M. & Lowe, R. Climate and non-climate drivers of dengue epidemics in Southern Coastal Ecuador. Am. J. Trop. Med. Hyg. 88, 971–981 (2013).
PubMed  PubMed Central  Article  Google Scholar 

7.
Jury, M. R. Climate influence on dengue epidemics in Puerto Rico. Int. J. Environ. Health Res. 18, 323–334 (2008).
PubMed  Article  PubMed Central  Google Scholar 

8.
Campbell, K. M. et al. Weather regulates location, timing, and intensity of dengue virus transmission between humans and mosquitoes. PLoS Negl. Trop. Dis. 9, e0003957 (2015).

9.
Adde, A. et al. Predicting dengue fever outbreaks in French Guiana using climate indicators. PLoS Negl. Trop. Dis. 10, e0004681 (2016).

10.
Dhimal, M. et al. Risk factors for the presence of chikungunya and dengue vectors (Aedes aegypti and Aedes albopictus), their altitudinal distribution and climatic determinants of their abundance in Central Nepal. PLoS Negl. Trop. Dis. 9, e0003545 (2015).
PubMed  PubMed Central  Article  CAS  Google Scholar 

11.
Descloux, E. et al. Climate-based models for understanding and forecasting dengue epidemics. PLoS Negl. Trop. Dis. 6, e1470 (2012).

12.
Aswi, A., Cramb, S. M., Moraga, P. & Mengersen, K. Epidemiology and infection Bayesian spatial and spatio-temporal approaches to modelling dengue fever: a systematic review. Epidemiol. Infect. 147, https://doi.org/10.1017/S0950268818002807 (2018).

13.
Johansson, M. A. et al. An open challenge to advance probabilistic forecasting for dengue epidemics. Proc. Natl Acad. Sci. USA 116, 24268–24274 (2019).

14.
Michael, E. et al. Continental-scale, data-driven predictive assessment of eliminating the vector-borne disease, lymphatic filariasis, in sub-Saharan Africa by 2020. BMC Med. 15, 176 (2017).
PubMed  PubMed Central  Article  Google Scholar 

15.
Smith, T. et al. Towards a comprehensive simulation model of malaria epidemiology and control. Parasitology 135, 1507–1516 (2008).

16.
Ryan, S. J., Carlson, C. J., Mordecai, E. A. & Johnson, L. R. Global expansion and redistribution of Aedes-borne virus transmission risk with climate change. PLoS Negl. Trop. Dis. 13, e0007213 (2019).

17.
Kraemer, M. U. et al. The global distribution of the arbovirus vectors Aedes aegypti and Ae. albopictus. eLife 4, e08347 (2015).

18.
Powell, J. R., Tabachnick, W. J., Powell, J. R. & Tabachnick, W. J. History of domestication and spread of Aedes aegypti—a review. Mem. Inst. Oswaldo Cruz. 108, 11–17 (2013).
PubMed  PubMed Central  Article  Google Scholar 

19.
Mordecai, E. A. et al. Detecting the impact of temperature on transmission of Zika, dengue, and chikungunya using mechanistic models. PLoS Negl. Trop. Dis. 11, e0005568 (2017).

20.
Shocket, M. S., Ryan, S. J. & Mordecai, E. A. Temperature explains broad patterns of Ross River virus transmission. eLife 7, e37762 (2018).

21.
Paull, S. H. et al. Drought and immunity determine the intensity of West Nile virus epidemics and climate change impacts. Proc. R. Soc. B Biol. Sci. 284, 20162078 (2017).
Article  Google Scholar 

22.
Costa EAP de, A., Santos EM de, M., Correia, J. C. & Albuquerque de, C. M. R. Impact of small variations in temperature and humidity on the reproductive activity and survival of Aedes aegypti (Diptera, Culicidae). Rev. Bras. Entomol. 54, 488–493 (2010).
Article  Google Scholar 

23.
Gaaboub, I. A., El-Sawaf, S. K. & El-Latif, M. A. Effect of different relative humidities and temperatures on egg-production and longevity of adults of Anopheles (Myzomyia) pharoensis Theob.1. Z. f.ür. Angew. Entomol. 67, 88–94 (2009).
Article  Google Scholar 

24.
Koenraadt, C. J. M. & Harrington, L. C. Flushing effect of rain on container-inhabiting mosquitoes Aedes aegypti and Culex pipiens (Diptera: Culicidae). J. Med. Entomol. 45, 28–35 (2009).
Google Scholar 

25.
Paaijmans, K. P., Wandago, M. O., Githeko, A. K., Takken, W. & Vulule, J. Unexpected high losses of Anopheles gambiae larvae due to rainfall. PLoS ONE 2, e1146 (2007).

26.
Benedum, C. M., Seidahmed, O. M. E., Eltahir, E. A. B. & Markuzon, N. Statistical modeling of the effect of rainfall flushing on dengue transmission in Singapore. PLoS Negl. Trop. Dis. 12, e0006935 (2018).

27.
Stewart Ibarra, A. M. et al. Dengue vector dynamics (Aedes aegypti) influenced by climate and social factors in Ecuador: implications for targeted control. PLoS ONE 8, e78263 (2013).

28.
Pontes, R. J., Spielman, A., Oliveira-Lima, J. W., Hodgson, J. C. & Freeman, J. Vector densities that potentiate dengue outbreaks in a Brazilian city. Am. J. Trop. Med Hyg. 62, 378–383 (2000).
CAS  PubMed  Article  PubMed Central  Google Scholar 

29.
Anyamba, A. et al. Climate teleconnections and recent patterns of human and animal disease outbreaks. PLoS Negl. Trop. Dis. 6, e1465 (2012).

30.
Huber, J. H., Childs, M. L., Caldwell, J. M. & Mordecai, E. A. Seasonal temperature variation influences climate suitability for dengue, chikungunya, and Zika transmission. PLoS Negl. Trop. Dis. 12, e0006451 (2018).

31.
Stewart-Ibarra, A. M. et al. Spatiotemporal clustering, climate periodicity, and social-ecological risk factors for dengue during an outbreak in Machala, Ecuador, in 2010. BMC Infect. Dis. 14, 610 (2014).
PubMed  PubMed Central  Article  Google Scholar 

32.
Agha, S. B., Tchouassi, D. P., Turell, M. J., Bastos, A. D. S. & Sang, R. Entomological assessment of dengue virus transmission risk in three urban areas of Kenya. PLoS Negl. Trop. Dis. 13, e0007686 (2019).

33.
Agha, S. B., Tchouassi, D. P., Bastos, A. D. S. & Sang, R. Dengue and yellow fever virus vectors: seasonal abundance, diversity and resting preferences in three Kenyan cities. Parasit. Vectors 10, 628 (2017).
PubMed  PubMed Central  Article  Google Scholar 

34.
Chretien, J.-P. et al. Drought-associated chikungunya emergence along coastal East Africa. Am. J. Trop. Med. Hyg. 76, 405–407 (2007).
PubMed  PubMed Central  Article  Google Scholar 

35.
Vu, D. M. et al. Unrecognized dengue virus infections in children, Western Kenya, 2014–2015. Emerg. Infect. Dis. 23, 1915–1917 (2017).
PubMed  PubMed Central  Article  Google Scholar 

36.
Gubler, D. J., Nalim, S., Saroso, J. S., Saipan, H. & Tan, R. Variation in susceptibility to oral infection with dengue viruses among geographic strains of Aedes Aegypti *. Am. J. Trop. Med. Hyg. 28, 1045–1052 (1979).
CAS  PubMed  Article  Google Scholar 

37.
Xavier-Carvalho, C., Chester Cardoso, C., de Souza Kehdya, F., Guilherme Pacheco, A. & Ozório Moraesa, M. Host genetics and dengue fever. Infect. Genet. Evol. 56, 99–110 (2017).
PubMed  Article  Google Scholar 

38.
Lourenço, J. & Recker, M. The 2012 Madeira Dengue Outbreak: epidemiological determinants and future epidemic potential. PLoS Negl. Trop. Dis. 8, e3083 (2014).

39.
Li, R. et al. Climate-driven variation in mosquito density predicts the spatiotemporal dynamics of dengue. Proc. Natl Acad. Sci. USA 119, 3624–3629 (2019).
Article  CAS  Google Scholar 

40.
Wang, X., Tang, S. & Cheke, R. A. A stage structured mosquito model incorporating effects of precipitation and daily temperature fluctuations. J. Theor. Biol. 411, 27–36 (2016).
PubMed  MATH  Article  PubMed Central  Google Scholar 

41.
Siraj, A. S. et al. Temperature modulates dengue virus epidemic growth rates through its effects on reproduction numbers and generation intervals. PLoS Negl. Trop. Dis. 11, e0005797 (2017).

42.
Oidtman, R. J. et al. Inter-annual variation in seasonal dengue epidemics driven by multiple interacting factors in Guangzhou, China. Nat. Commun. 10, 1–12 (2019).

43.
Pyper, B. J. & Peterman, R. M. Comparison of methods to account for autocorrelation in correlation analyses of fish data. Can. J. Fish. Aquat. Sci. 55, 2127–2140 (1998).
Article  Google Scholar 

44.
Hurtado-Daz, M., Riojas-Rodrguez, H., Rothenberg, S., Gomez-Dantes, H. & Cifuentes, E. Impact of climate variability on the incidence of dengue in Mexico. Trop. Med. Int. Heal. 12, 1327–1337 (2007).

45.
Mordecai, E. A. et al. Thermal biology of mosquito‐borne disease. Ecol. Lett. 22, 1690–1708 (2019).

46.
Carrington, L. B., Armijos, M. V., Lambrechts, L., Barker, C. M. & Scott, T. W. Effects of fluctuating daily temperatures at critical thermal extremes on Aedes aegypti life-history traits. PLoS ONE 8, e58824 (2013).

47.
Ngugi, H. N. et al. Characterization and productivity profiles of Aedes aegypti (L.) breeding habitats across rural and urban landscapes in western and coastal Kenya. Parasit. Vectors 10, 331 (2017).
PubMed  PubMed Central  Article  Google Scholar 

48.
Lowe, R. et al. Nonlinear and delayed impacts of climate on dengue risk in Barbados: a modelling study. PLoS Med. 15, e1002613 (2018).

49.
Laureano-Rosario, A. E., Garcia-Rejon, J. E., Gomez-Carro, S., Farfan-Ale, J. A. & Muller-Kargera, F. E. Modelling dengue fever risk in the State of Yucatan, Mexico using regional-scale satellite-derived sea surface temperature. Acta Trop. 172, 50–57 (2017).
PubMed  Article  PubMed Central  Google Scholar 

50.
Li, C. et al. Modeling and projection of dengue fever cases in Guangzhou based on variation of weather factors. Sci. Total Environ. 605–606, 867–873 (2017).
ADS  PubMed  Article  CAS  PubMed Central  Google Scholar 

51.
Li, C. F., Lim, T. W., Han, L. L. & Fang, R. Rainfall, abundance of Aedes aegypti and dengue infection in Selangor, Malaysia. Southeast Asian J. Trop. Med Public Health 16, 560–568 (1985).
CAS  PubMed  PubMed Central  Google Scholar 

52.
Johansson, M. A., Dominici, F. & Glass, G. E. Local and global effects of climate on dengue transmission in Puerto Rico. PLoS Negl. Trop. Dis. 3, e382 (2009).

53.
Kenneson, A. et al. Social-ecological factors and preventive actions decrease the risk of dengue infection at the household-level: Results from a prospective dengue surveillance study in Machala, Ecuador. PLoS Negl. Trop. Dis. 11, e0006150 (2017).

54.
Reich, N. G. et al. Interactions between serotypes of dengue highlight epidemiological impact of cross-immunity. J. R. Soc. Interface 10, 20130414 (2013).
PubMed  PubMed Central  Article  Google Scholar 

55.
Wen, J. et al. Dengue virus-reactive CD8+ T cells mediate cross-protection against subsequent Zika virus challenge. Nat. Commun. 8, 1459 (2017).
ADS  PubMed  PubMed Central  Article  CAS  Google Scholar 

56.
Rodriguez-Barraquer, I., Salje, H. & Cummings, D. A. Opportunities for improved surveillance and control of dengue from age-specific case data. eLife 8, e45474 (2019).

57.
Stoddard, S. T. et al. House-to-house human movement drives dengue virus transmission. Proc. Natl Acad. Sci. USA 110, 994–999 (2013).
ADS  CAS  PubMed  Article  Google Scholar 

58.
Wesolowski, A. et al. Impact of human mobility on the emergence of dengue epidemics in Pakistan. Proc. Natl Acad. Sci. USA 112, 11887–11892 (2015).
ADS  CAS  PubMed  Article  Google Scholar 

59.
Vaidya, A., Bravo-Salgado, A. D. & Mikler, A. R.. Modeling climate-dependent population dynamics of mosquitoes to guide public health policies. in Proceedings of the 5th ACM Conference on Bioinformatics, Computational Biology, and Health Informatics Vol. 14 (eds Baldi, P. & Wang, W.) 380–389 (Newport Beach, CA, USA, 2014).

60.
Schmidt, C. A., Comeau, G., Monaghan, A. J., Williamson, D. J. & Ernst, K. C. Effects of desiccation stress on adult female longevity in Aedes aegypti and Ae. albopictus (Diptera: Culicidae): results of a systematic review and pooled survival analysis. Parasit. Vectors 11, 267 (2018).
PubMed  PubMed Central  Article  Google Scholar 

61.
Vazquez-Prokopec, G. M., Galvin, W. A., Kelly, R. & Kitron, U. A new, cost-effective, battery-powered aspirator for adult mosquito collections. J. Med. Entomol. 46, 1256–1259 (2009).
PubMed  PubMed Central  Article  Google Scholar 

62.
Waggoner, J. J. et al. Single-reaction multiplex reverse transcription PCR for detection of Zika, chikungunya, and dengue viruses. Emerg. Infect. Dis. 22, 1295–1297 (2016).
CAS  PubMed  PubMed Central  Article  Google Scholar 

63.
Lanciotti, R. S., Calisher, C. H., Gubler, D. J., Chang, G. J. & Vorndam, A. V. Rapid detection and typing of dengue viruses from clinical samples by using reverse transcriptase-polymerase chain reaction. J. Clin. Microbiol. 30, 545–551 (1992).
CAS  PubMed  PubMed Central  Article  Google Scholar 

64.
Grossi-Soyster, E. N. et al. Serological and spatial analysis of alphavirus and flavivirus prevalence and risk factors in a rural community in western Kenya. PLoS Negl. Trop. Dis. 11, e0005998 (2017).

65.
Palamara, G. M. et al. Inferring the temperature dependence of population parameters: the effects of experimental design and inference algorithm. Ecol. Evol. 4, 4736–4750 (2014).
PubMed  PubMed Central  Article  Google Scholar 

66.
Team R. C. R.: A language and environment for statistical computing. R Found Stat. Comput. https://www.r-project.org (2018).

67.
Shocket, M. S. et al. Environmental drivers of vector-borne disease. in Population Biology of Vector-borne Diseases. (eds Drake, J. M., Bonsall, M. B. & Strand, M. R.) 85–118 (Oxford University Press, 2020).

68.
Colón-González, F. J., Bentham, G. & Lake, I. R. Climate variability and dengue fever in warm and humid Mexico. Am. J. Trop. Med. Hyg. 84, 757–763 (2011).
PubMed  PubMed Central  Article  Google Scholar 

69.
Wang, C., Jiang, B., Fan, J., Wang, F. & Liu, Q. A study of the dengue epidemic and meteorological factors in Guangzhou, China, by using a zero-inflated Poisson regression model. Asia Pac. J. Public Heal. 26, 48–57 (2014).
Article  Google Scholar 

70.
Minh An, D. T. & Rocklöv, J. Epidemiology of dengue fever in Hanoi from 2002 to 2010 and its meteorological determinants. Glob. Health Action. 7, 23074 (2014).
PubMed  Article  PubMed Central  Google Scholar 

71.
Wu, P.-C., Guoa, H.-R., Lung, S.-C., Lin, C.-Y. & Su, H.-J. Weather as an effective predictor for occurrence of dengue fever in Taiwan. Acta Trop. 103, 50–57 (2007).
PubMed  Article  PubMed Central  Google Scholar 

72.
Karim, M. N., Munshi, S. U., Anwar, N. & Alam, M. S. Climatic factors influencing dengue cases in Dhaka city: a model for dengue prediction. Indian J. Med. Res. 136, 32–39 (2012).
PubMed  PubMed Central  Google Scholar 

73.
Nakhapakorn, K. & Tripathi, N. An information value based analysis of physical and climatic factors affecting dengue fever and dengue haemorrhagic fever incidence. Int. J. Health Geogr. 4, 13 (2005).
PubMed  PubMed Central  Article  Google Scholar 

74.
Gharbi, M. et al. Time series analysis of dengue incidence in Guadeloupe, French West Indies: forecasting models using climate variables as predictors. BMC Infect. Dis. 11, 166 (2011).
PubMed  PubMed Central  Article  Google Scholar 

75.
Sharmin, S., Glass, K., Viennet, E. & Harley, D. Interaction of mean temperature and daily fluctuation influences dengue incidence in Dhaka, Bangladesh. PLoS Negl. Trop. Dis. 9, e0003901 (2015).

76.
Sriprom, M., Chalvet-Monfray, K., Chaimane, T., Vongsawat, K. & Bicout, D. J. Monthly district level risk of dengue occurrences in Sakon Nakhon Province, Thailand. Sci. Total Environ. 408, 5521–5528 (2010).
ADS  CAS  PubMed  Article  PubMed Central  Google Scholar 

77.
Martínez-Bello, D., López-Quílez, A. & Prieto, A. T. Spatiotemporal modeling of relative risk of dengue disease in Colombia. Stoch. Environ. Res. Risk Assess. 32, 1587–1601 (2018).
Article  Google Scholar 

78.
Didan, K., Barreto Munoz, A., Solano, R. & Huete, A. MODIS vegetation index user’s guide (MOD13 Series) [Internet]. https://vip.arizona.edu/documents/MODIS/MODIS_VI_UsersGuide_June_2015_C6.pdf (2015).

79.
Sulla-Menashe, D. & Friedl, M. A. User guide to collection 6 MODIS land cover (MCD12Q1 and MCD12C1) product. https://icdc.cen.uni-hamburg.de/fileadmin/user_upload/icdc_Dokumente/MODIS/mcd12_user_guide_v6.pdf (2018). More

1.
Barbier, E. B. et al. The value of estuarine and coastal ecosystem services. Ecol. Monogr. 81, 169–193 (2011).
Article  Google Scholar 
2.
Costanza, R. et al. Changes in the global value of ecosystem services.Glob. Environ. Chang. 26, 152–158 (2014).
Article  Google Scholar 

3.
Airoldi, L. & Beck, M. W. Loss, status and trends for coastal marine habitats of Europe. Oceanogr. Mar. Biol. Annu. Rev. 45, 345–405 (2007).
Google Scholar 

4.
Kainuma, Mami et al. Current status of mangroves worldwide. Middle East 624, 0–4 (2013).
Google Scholar 

5.
Fagherazzi, S. et al. Sea level rise and the dynamics of the marsh-upland boundary. Front. Environ. Sci. 7, 25 (2019).
Article  Google Scholar 

6.
Kirwan, M. L. & Gedan, K. B. Sea-level driven land conversion and the formation of ghost forests. Nat. Clim. Change 9, 450–457 (2019).
Article  Google Scholar 

7.
Craft, C. et al. Forecasting the effects of accelerated sea‐level rise on tidal marsh ecosystem services. Front. Ecol. Environ. 7, 73–78 (2009).
Article  Google Scholar 

8.
Nicholls, R. J. & Cazenave, A. Sea-level rise and its impact on coastal zones. Science 328, 1517–1520 (2010).
CAS  Article  Google Scholar 

9.
Schuerch, M. et al. Modeling the influence of changing storm patterns on the ability of a salt marsh to keep pace with sea level rise. J. Geophys. Res. Earth Surf. 118, 84–96 (2013).
Article  Google Scholar 

10.
Temmerman, S. et al. Ecosystem-based coastal defence in the face of global change. Nature 504, 79–83 (2013).
CAS  Article  Google Scholar 

11.
Syvitski, J. P. et al. Impact of humans on the flux of terrestrial sediment to the global coastal ocean. Science 308, 376–380 (2005).
CAS  Article  Google Scholar 

12.
Ezcurra, E. et al. A natural experiment reveals the impact of hydroelectric dams on the estuaries of tropical rivers.Sci. Adv. 5, eaau9875 (2019).
CAS  Article  Google Scholar 

13.
Kirwan, M. L. et al. Overestimation of marsh vulnerability to sea level rise. Nat. Clim. Change 6, 253–260 (2016).
Article  Google Scholar 

14.
Schuerch, M. et al. Future response of global coastal wetlands to sea-level rise. Nature 561, 231–234 (2018).
CAS  Article  Google Scholar 

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

16.
Gittman, R. K., Scyphers, S. B., Smith, C. S., Neylan, I. P. & Grabowski, J. H. Ecological consequences of shoreline hardening: a meta-analysis. BioScience 66, 763–773 (2016).
Article  Google Scholar 

17.
Smith, C. S. et al. Hurricane damage along natural and hardened estuarine shorelines: Using homeowner experiences to promote nature-based coastal protection. Mar. Policy 81, 350–358 (2017).
Article  Google Scholar 

18.
Shepard, C. C., Crain, C. M. & Beck, M. W. The protective role of coastal marshes: a systematic review and meta-analysis. PLoS ONE 6, e27374 (2011).
CAS  Article  Google Scholar 

19.
Gedan, K. B., Kirwan, M. L., Wolanski, E., Barbier, E. B. & Silliman, B. R. The present and future role of coastal wetland vegetation in protecting shorelines: answering recent challenges to the paradigm. Clim. Change 106, 7–29 (2011).
Article  Google Scholar 

20.
Leonardi, N., Ganju, N. K. & Fagherazzi, S. A linear relationship between wave power and erosion determines salt-marsh resilience to violent storms and hurricanes. Proc. Nat. Acad. Sci. USA 113, 64–68 (2016).
CAS  Article  Google Scholar 

21.
Barbier, E. B. et al. Coastal ecosystem-based management with nonlinear ecological functions and values. Science 319, 321–323 (2008).
CAS  Article  Google Scholar 

22.
Cohen-Shacham, E., Walters, G., Janzen, C. & Maginnis, S. Nature-based Solutions to Address Global Societal Challenges (IUCN, 2016).

23.
Fargione, J. E. et al. Natural climate solutions for the United States. Sci. Adv. 4, eaat1869 (2018).
Article  Google Scholar 

24.
Seddon, N. et al. Global recognition of the importance of Nature-based Solutions to the impacts of climate change. Glob. Sustain. 3, 1–12 (2020).
Article  Google Scholar 

25.
Bilkovic, D. M. et al. Living Shorelines: The Science and Management of Nature-Based Coastal Protection (CRC Press, 2017).

26.
Bayraktarov, E. et al. The cost and feasibility of marine coastal restoration. Ecol. Appl. 26, 1055–1074 (2016).
Article  Google Scholar 

27.
Liu, Z., Cui, B. & He, Q. Shifting paradigms in coastal restoration: Six decades’ lessons from China. Sci. Total Environ. 566, 205–214 (2016).
Article  CAS  Google Scholar 

28.
Turner, R. K., Burgess, D., Hadley, D., Coombes, E. & Jackson, N. A cost–benefit appraisal of coastal managed realignment policy.Glob. Environ. Chang. 17, 397–407 (2007).
Article  Google Scholar 

29.
Donatelli, C., Ganju, N. K., Zhang, X., Fagherazzi, S. & Leonardi, N. Salt marsh loss affects tides and the sediment budget in shallow bays. J. Geophys. Res. Earth Surf. 123, 2647–2662 (2018).
Article  Google Scholar 

30.
Benayas, J. M. R., Newton, A. C., Diaz, A. & Bullock, J. M. Enhancement of biodiversity and ecosystem services by ecological restoration: a meta-analysis. Science 325, 1121–1124 (2009).
CAS  Article  Google Scholar 

31.
Friess, D. A. et al. Are all intertidal wetlands naturally created equal? Bottlenecks, thresholds and knowledge gaps to mangrove and saltmarsh ecosystems. Biol. Rev. 87, 346–366 (2012).
Article  Google Scholar 

32.
Webb, E. L. et al. A global standard for monitoring coastal wetland vulnerability to accelerated sea-level rise. Nature Clim. Change 3, 458–465 (2013).
Article  Google Scholar 

33.
Hu, Z. et al. Revegetation of a native species in a newly formed tidal marsh under varying hydrological conditions and planting densities in the Yangtze Estuary. Ecol. Eng. 83, 354–363 (2015).
Article  Google Scholar 

34.
Phillips, D. H. et al. Impacts of mangrove density on surface sediment accretion, belowground biomass and biogeochemistry in Puttalam Lagoon, Sri Lanka. Wetlands 37, 471–483 (2017).
Article  Google Scholar 

35.
Kirwan, M. L. et al. Limits on the adaptability of coastal marshes to rising sea level. Geophys. Res. Lett. 37, L23401 (2010).
Article  Google Scholar 

36.
Turner, R. E., Baustian, J. J., Swenson, E. M. & Spicer, J. S. Wetland sedimentation from hurricanes Katrina and Rita. Science 314, 449–452 (2006).
CAS  Article  Google Scholar 

37.
French, C. E., French, J. R., Clifford, N. J. & Watson, C. J. Sedimentation-erosion dynamics of abandoned reclamations: the role of waves and tides. Cont. Shelf Res. 20, 1711–1733 (2000).
Article  Google Scholar 

38.
Cahoon, D. R. et al. High-precision measurements of wetland sediment elevation: II. The rod surface elevation table. J. Sediment. Res. 72, 734–739 (2002).
CAS  Article  Google Scholar 

39.
Cahoon, D. R. A review of major storm impacts on coastal wetland elevations. Estuar. Coast. 29, 889–898 (2006).
Article  Google Scholar 

40.
Howe, A. J., Rodriguez, J. F. & Saco, P. M. Surface evolution and carbon sequestration in disturbed and undisturbed wetland soils of the Hunter estuary, southeast Australia. Estuar. Coast. Shelf Sci. 84, 75–83 (2009).
CAS  Article  Google Scholar 

41.
Krauss, K. W. et al. Created mangrove wetlands store belowground carbon and surface elevation change enables them to adjust to sea-level rise. Sci. Rep. 7, 1–11 (2017).
Article  CAS  Google Scholar 

42.
Carey, J. C., Moran, S. B., Kelly, R. P., Kolker, A. S. & Fulweiler, R. W. The declining role of organic matter in New England salt marshes. Estuar. Coast 40, 626–639 (2017).
CAS  Article  Google Scholar 

43.
Lovelock, C. E. et al. The vulnerability of Indo-Pacific mangrove forests to sea-level rise. Nature 526, 559–563 (2015).
CAS  Article  Google Scholar 

44.
Anisfeld, S. C., Hill, T. D. & Cahoon, D. R. Elevation dynamics in a restored versus a submerging salt marsh in Long Island Sound. Estuar. Coast. Shelf Sci. 170, 145–154 (2016).
Article  Google Scholar 

45.
Baustian, J. J., Mendelssohn, I. A. & Hester, M. W. Vegetation’s importance in regulating surface elevation in a coastal salt marsh facing elevated rates of sea level rise. Glob. Chang. Biol. 18, 3377–3382 (2012).
Article  Google Scholar 

46.
Cahoon, D. R., French, J. R., Spencer, T., Reed, D. & Möller, I. Vertical accretion versus elevational adjustment in UK saltmarshes: an evaluation of alternative methodologies. Geol. Soc. Lond. Spec. Publ. 175, 223–238 (2000).
Article  Google Scholar 

47.
Spencer, T. et al. Surface elevation change in natural and re-created intertidal habitats, eastern England, UK, with particular reference to Freiston Shore. Wetl. Ecol. Manag. 20, 9–33 (2012).
Article  Google Scholar 

48.
Craft, C. et al. The pace of ecosystem development of constructed Spartina alterniflora marshes. Ecol. Appl. 13, 1417–1432 (2003).
Article  Google Scholar 

49.
Duarte, C. M., Losada, I. J., Hendriks, I. E., Mazarrasa, I. & Marbà, N. The role of coastal plant communities for climate change mitigation and adaptation. Nat. Clim. Change 3, 961–968 (2013).
CAS  Article  Google Scholar 

50.
Fagherazzi, S. et al. Numerical models of salt marsh evolution: Ecological, geomorphic, and climatic factors. Rev. Geophys. 50, RG1002 (2012).
Article  Google Scholar 

51.
Smith, C. S., Puckett, B., Gittman, R. K. & Peterson, C. H. Living shorelines enhanced the resilience of saltmarshes to Hurricane Matthew. Ecol. Appl. 28, 871–877 (2018).
Article  Google Scholar 

52.
Oosterlee, L. et al. Tidal marsh restoration design affects feedbacks between inundation and elevation change. Estuar. Coast. 41, 613–625 (2018).
Article  Google Scholar 

53.
Ganju, N. K. Marshes are the new beaches: integrating sediment transport into restoration planning. Estuar. Coast. 42, 917–926 (2019).
CAS  Article  Google Scholar 

54.
Ford, M. A., Cahoon, D. R. & Lynch, J. C. Restoring marsh elevation in a rapidly subsiding salt marsh by thin-layer deposition of dredged material. Ecol. Eng. 12, 189–205 (1999).
Article  Google Scholar 

55.
Temmerman, S., Govers, G., Wartel, S. & Meire, P. Spatial and temporal factors controlling short‐term sedimentation in a salt and freshwater tidal marsh, Scheldt estuary, Belgium, SW Netherlands. Earth Surf. Processes Landforms 28, 739–755 (2003).
Article  Google Scholar 

56.
Morris, J. T., Sundareshwar, P. V., Nietch, C. T., Kjerfve, B. & Cahoon, D. R. Responses of coastal wetlands to rising sea level. Ecology 83, 2869–2877 (2002).
Article  Google Scholar 

57.
Mudd, S. M., D’Alpaos, A. & Morris, J. T. How does vegetation affect sedimentation on tidal marshes? Investigating particle capture and hydrodynamic controls on biologically mediated sedimentation. J. Geophys. Res. Earth Surf. 115, F03029 (2010).
Google Scholar 

58.
Fricke, A. T., Nittrouer, C. A., Ogston, A. S. & Vo-Luong, H. P. Asymmetric progradation of a coastal mangrove forest controlled by combined fluvial and marine influence, Cù Lao Dung, Vietnam. Cont. Shelf Res. 147, 78–90 (2017).
Article  Google Scholar 

59.
Möller, I., Spencer, T., French, J. R., Leggett, D. J. & Dixon, M. Wave transformation over salt marshes: a field and numerical modelling study from North Norfolk, England. Estuar. Coast. Shelf Sci. 49, 411–426 (1999).
Article  Google Scholar 

60.
Jadhav, R. S., Chen, Q. & Smith, J. M. Spectral distribution of wave energy dissipation by salt marsh vegetation. Coast. Eng. 77, 99–107 (2013).
Article  Google Scholar 

61.
Kirwan, M. L. & Guntenspergen, G. R. Influence of tidal range on the stability of coastal marshland. J. Geophys. Res. Earth Surf. 115, F02009 (2010).
Article  Google Scholar 

62.
Ganju, N. K., Nidzieko, N. J. & Kirwan, M. L. Inferring tidal wetland stability from channel sediment fluxes: Observations and a conceptual model. J. Geophys. Res. Earth Surf. 118, 2045–2058 (2013).
Article  Google Scholar 

63.
Zhang, X. et al. Determining the drivers of suspended sediment dynamics in tidal marsh-influenced estuaries using high-resolution ocean color remote sensing. Remote Sens. Environ. 240, 111682 (2020).
Article  Google Scholar 

64.
Hopkinson, C. S., Morris, J. T., Fagherazzi, S., Wollheim, W. M. & Raymond, P. A. Lateral marsh edge erosion as a source of sediments for vertical marsh accretion. J. Geophys. Res. Biogeo. 123, 2444–2465 (2018).
CAS  Article  Google Scholar 

65.
Castagno, K. A. et al. Intense storms increase the stability of tidal bays. Geophys. Res. Lett. 45, 5491–5500 (2018).
Article  Google Scholar 

66.
Walling, D. E. The Impact of Global Change on Erosion and Sediment Transport by Rivers: Current Progress and Future Challenges (UNESCO, 2009).

67.
Yu, Y. et al. New discharge regime of the Huanghe (Yellow River): causes and implications. Cont. Shelf Res. 69, 62–72 (2013).
Article  Google Scholar 

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

69.
Donatelli, C., Kalra, T. S., Fagherazzi, S., Zhang, X. & Leonardi, N. Dynamics of marsh‐derived sediments in lagoon‐type estuaries. J. Geophys. Res. Earth Surf. 125, e2020JF005751 (2020).
Article  Google Scholar 

70.
Peteet, D. M. et al. Sediment starvation destroys New York City marshes’ resistance to sea level rise. Proc. Nat. Acad. Sci. USA 115, 10281–10286 (2018).
CAS  Article  Google Scholar 

71.
Reed, D. J. Understanding tidal marsh sedimentation in the Sacramento-San Joaquin Delta, California. J. Coastal Res. 36, 605–611 (2002).
Article  Google Scholar 

72.
Cahoon, D. R., Lynch, J. C., Roman, C. T., Schmit, J. P. & Skidds, D. E. Evaluating the relationship among wetland vertical development, elevation capital, sea-level rise, and tidal marsh sustainability. Estuar. Coast. 42, 1–15 (2019).
CAS  Article  Google Scholar 

73.
Kondolf, G. M., Rubin, Z. K. & Minear, J. T. Dams on the Mekong: Cumulative sediment starvation. Water Resour. Res. 50, 5158–5169 (2014).
Article  Google Scholar 

74.
Reed, D. J. Patterns of sediment deposition in subsiding coastal salt marshes, Terrebonne Bay, Louisiana: the role of winter storms. Estuaries 12, 222–227 (1989).
Article  Google Scholar 

75.
Ganju, N. K. et al. Spatially integrative metrics reveal hidden vulnerability of microtidal salt marshes. Nat. Commun. 8, 14156 (2017).
CAS  Article  Google Scholar 

76.
Vörösmarty, C. J. et al. Anthropogenic sediment retention: major global impact from registered river impoundments. Glob. Planet. Change 39, 169–190 (2003).
Article  Google Scholar 

77.
Borenstein, M., Hedges, L. V., Higgins, J. P. T. & Rothstein, H. R. Introduction to Meta-Analysis (John Wiley & Sons, Ltd., 2009). More

Experimental protocol
Twenty healthy male adults participated in the study. To limit the effects of differences in the participants’ anthropometry, the volunteers were matched according to their height and weight.
The leg and hand dominances of each participant was noted before running the experiment. The subjects were placed at random on the left ((#1—right hand holding the load) or the right (#2—left hand holding the load) side of the carried object. It happened by chance that three lefthanders were affected to #2. The individuals had an average height of (mean ± SD) 1.77 ± 0.07 m (#1), and 1.77 ± 0.05 m (#2), and an average weight of 74.78 ± 9.00 kg (#1), and 74.54 ± 12.38 kg (#2). The load was symmetrical in shape, and its weight was evenly balanced between the participants who were positioned randomly with respect to the object transported to counteract the effect of a dominant side (Fig. 1A).
The study was carried out with healthy individuals who wrote their informed consent to participate in the experiment and to be filmed and photographed. The experiment was non-interventional, and the movements performed by the volunteers were no more risky than those they perform in daily activities. The study was approved by the Research Ethics Committee of the University of Toulouse, France (Number IRB00011835-2019-11-26-172, Université Fédérale de Toulouse IRB #1).
The instructions given to the volunteers were: “Move the load together from point A to point B” and “Any communication between you is forbidden during the experiment”. Point A and point B laid 20 m apart. No explicit instruction was given as to how fast the volunteers should perform the task. The volunteers were tested with three conditions called CT20, CT30, and CT40 corresponding to a load representing on average 20%, 30%, and 40% of the sum of their body masses respectively. The three conditions were tested in random order for each pair. To avoid adaptations due to familiarization or learning, only one trial per pairs and condition was recorded.
Kinematics and kinetics
Thirteen MX3, and TS40 Vicon cameras (Vicon©, Oxford) were used to capture the positions of ninety-one retro-reflective markers taped on the system formed by the paired individuals and the load they carry (hereafter called Poly-Articulated Collective System—PACS): 42 markers on each individual31,32, and seven on the load (Fig. 1A). The acquisition frequency was set to 200 Hz. In order to record the walking patterns of the individuals at a stable speed, and thus to exclude the acceleration and deceleration phases at the beginning and end of each trial, the calibrated volume (30m3) corresponded to the central part of the walkway. This covered about two steps. Concerning kinematic analysis, the PACS was reconstructed with the Vicon Nexus™ 1.8.5 software. Reconstruction was impossible for one pair of individuals who had lost one reflective marker. The two lateral handles on each side of the load were equipped with a 6-axis force sensor (Sensix®, France) (Fig. 1B), allowing to record the reaction forces and moments at a sampling frequency of 2000 Hz. The kinematic and kinetic measurement errors were 1 mm for 1 m for the positions (Vicon system) and ± 0.01 N for the forces (Sensix sensors), respectively. The sensors frames were located with the help of screwed reflective markers. The data were filtered with 4th order Butterworth filters with a cut-off frequency of 5 Hz for the kinematic data, and of 10 Hz for the kinetic data. To ensure at least one complete walking cycle for each subject of a pair, the gait cycle of the PACS was defined from the first heel strike of individual #1 to the third heel strike of individual #2.
COM determination and related parameters
The carried object, which constituted the 33th segment of the PACS, was built in aluminum and was therefore extremely rigid. It was completely symmetrical about its sagittal plane (Fig. 1B) and therefore its weight was evenly balanced between the participants. The object was equipped with a rod at its center where standard cast iron discs could be slid to increase its weight. The CoM of the object was determined at the intersection point of the vertical lines obtained by hanging the object without discs with a thread fixed at different positions. When the object was loaded, the position of its CoM was then adjusted by taking into account the added cast iron discs and by considering a homogeneous mass distribution inside the discs.
The De Leva Anthropometric table33 allowed us to estimate the mass mi as well as the CoM of each segment i (CoMi) of the PACS, and thus to compute its global CoM (CoMPACS) as follows:

$${varvec{G}}{text{PACS}} = frac{1}{{m_{PACS} }}mathop sum limits_{i = 1}^{n = 33} m_{i} {varvec{G}}_{i}$$
(1)

with GPACS corresponding to the 3D position of the CoMPACS in the frame R (the global coordinate system), mPACS to the mass of the PACS, n the number of PACS segments (i.e. 16 segments per volunteer, plus one segment for the box), and Gi corresponding to the 3D position of the CoMi in R.
The vertical amplitude (Az = Zmax − Zmin, in meters) of the CoMPACS trajectory along two consecutive steps, and the length of two consecutive steps by each individual were also computed.
Assessment of energetic exchanges
To assess energetic exchanges, forward kinetic work, vertical work, and external work of the forces applied to the CoMPACS were computed25.
Forward kinetic work (Wkf) was defined as the positive work needed to move the CoMPACS forward, and it was calculated as the sum of the increments of forward kinetic energy (Ekf) along the time curve:

$$E_{{{text{kf}}}} = { }frac{1}{2} m overrightarrow{V_{f}}left( {text{t}} right)^{{2}} _{/R}$$
(2)

with m being the mass of the individual, and (overrightarrow {{V_{f} }}) (t)/R the linear forward velocity of the CoMPACS in the frame R. The x-, y- and z-axis of the frame R, corresponding to the medio-lateral, antero-posterior, and vertical directions respectively, are illustrated in Fig. 1A.
Vertical work (Wv) was defined as the positive work needed to move the CoMPACS against gravity, and it was calculated as the sum of the increments of the vertical kinetic energy (Ekv) plus the potential energy (Epot) along the time curve with:

$$E_{{{text{kv}}}} = { }frac{1}{2} m overrightarrow{V_{v}}left( {text{t}} right)^{{2}} _{/R}$$
(3)

and

$$E_{{{text{pot}}}} = mgh_{{/{text{R}}}}$$
(4)

where (overrightarrow {{V_{v} }}) (t)/R is the linear vertical velocity of the CoMPACS in R, g = 9.81 m s−2 is the acceleration due to gravity, and h/R is the height of the CoMPACS in R.
The external work (Wext), corresponding to the positive external work needed to raise and accelerate the CoMPACS, was computed as the sum of the increments of the external mechanical energy (Eext) along the time curve with:

$$E_{{{text{ext}}}} = E_{{{text{pot}}}} + E_{{{text{kv}}}} + E_{{{text{kf}}}}$$
(5)

The energy recovered (called recovery rate (RR)10) by the CoMPACS in the sagittal plane was computed with the following formula17:

$$RR = { 1}00frac{{W{text{kf}} + W{text{v}} – W{text{ext}}}}{{W{text{kf}} + W{text{v}}}}$$
(6)

RR is the percentage of kinetic energy converted into potential energy7,24,25,34,35 and vice versa.
In the present study, internal work was also considered in order to encompass the coordination between all body segments. Based on the assumption of a conservative Poly-Articulated Model (PAM), internal work (Wint) was computed as the sum of the increments of the Eint,k along the time curve with:

$$E_{{{text{int}},{text{k}}}} = frac{1}{2}~mathop sum limits_{{i = 1}}^{{33}} (m_{i} overrightarrow {{V_{{~i}} }} left( {text{t}} right)^{{text{2}}} _{{/{text{R}}*}} + m_{{text{i}}} K_{{i^{2} }} {text{ }} times vec{omega }^{2} _{i} /_{{{text{R}}*}} )$$
(7)

where mi is the mass of the ith segment, (overrightarrow {{V_{i} }})(t)/R* the linear velocity of its CoM in the sagittal plane of the barycentric coordinate system (R*), Ki its radius of gyration around its CoM, and (vec{omega }_{i})2/R* its angular velocity in R* 36.
The total mechanical energy of the PACS (Etot) was computed as the sum of the internal kinetic energy (Eint,k) of each segment, plus the potential energy (Epot), and the forward (Ekf ) and vertical (Ekv ) kinetic energy of the CoMPACS in the sagittal plane21,25,37,38:

$$E_{{{text{tot}}}} = E_{{{text{int}},{text{k}}}} + E_{{{text{pot}}}} + E_{{{text{kf}}}} + E_{{{text{kv}}}}$$
(8)

Finally, the total mechanical power (PmecaTot) was used to assess the amount of energy spent or gained by the CoMPACS per unit of time (Δt):

$$P_{{{text{mecaTot}}}} = frac{{W{text{ext}}}}{Delta t} + frac{{{ }W{text{int}}}}{Delta t} = P_{{{text{ext}}}} + P_{{{text{int}}}}$$
(9)

Calculation of internal efforts
The resultant joint moments at the wrist, elbow, shoulder, neck, and back joints were calculated using a bottom-up Newton–Euler recursive algorithm39. Cardanic angles were used to represent the rotation of the segments coordinate system relative to the global coordinate system40. The segment masses, inertia tensors, and center of mass locations were estimated for each subject according to the scaling equations proposed in Dumas et al. (2007)41. In order to estimate the muscular torque produced at all the joints of the upper-limbs, shoulders, neck, and back, the Moment Cost Function (MCF in kg m2 s−2, 42) was computed as follows:

$${text{MCF}} = sqrt {M_{L_wt}^{2} } + sqrt {M_{R_wt}^{2} } + sqrt {M_{L_el}^{2} } + sqrt {M_{R_el}^{2} } + sqrt {M_{L_sh}^{2} } + sqrt {M_{R_sh}^{2} } + sqrt {M_{back}^{2} } + sqrt {M_{neck}^{2} }$$
(10)

where ML_wt, MR_wt, ML_el, MR_el, ML_sh, MR_sh, Mback, and Mneck are the mean values over a PACS gait cycle of the three-dimensional left and right wrist, left and right elbow, left and right shoulder, top of the back and neck moments, respectively. (sqrt {{text{M}}^{2} }) represents the Euclidian norm of M, i.e. (sqrt {sumnolimits_{i = 1}^{3} {left( {M_{i} } right)^{2} } }), with Mi the i-th component of the vector M.
We summed the MCF values of the two individuals of each pair to obtain the total moment cost function (TotMCF). The TotMCF allows to quantify the global muscular effort developed at the upper-limbs of the PACS during one gait cycle of the carrying of the load. The MCF difference (∆MCF) between the two individuals was also computed to investigate whether the volunteers developed the same efforts while carrying the object.
Kinetic synergy analysis
We extracted the synergies by using a principal component analysis (PCA) applied to the wrist, elbow, shoulder, back, and neck joint moment on the right and left sides of the body. The PCA was used to reduce data dimensionality13,35,43. It consisted in the eigen-decomposition of the co-variance matrix of the joint moment data (Matlab eig function). The joint moments data were arranged in time × joint moment matrices. We called the eigenvectors extracted from the PCA synergy vectors13. The number of synergies was determined from the VAF (Variance Accounted For), which corresponds to the cumulative sum of the eigenvalues, ordered from the greatest to the lowest value, normalized by the total variance computed as the sum of all eigenvalues. We defined the number of synergies as the first number for which the VAF was greater than 0.9. The synergy vectors retained were then rotated using a Varimax rotation method to improve interpretability44.
We extracted the synergy vectors for each experimental condition separately. We first performed an analysis on each individual separately. In this analysis the initial data matrices were constituted of all available time frames in line, concatenated, and of eight columns corresponding to each joint moment, namely the right wrist, left wrist, right elbow, left elbow, right shoulder, left shoulder, back, and neck. The values in the matrix corresponded to the norm of the joint moment vector at a given time frame. We then performed a second analysis to identify possible co-variations between the joint moments of the two participants in each pair. The columns of the initial matrices were thus constituted of the joint moments of the two loaded arms, i.e., the right wrist, elbow, and shoulder joint moments of participant #1, plus the left wrist, elbow and shoulder joint moments of participant #2. The synergy vectors were compared across conditions by computing Pearson’s r correlations on their PCA weightings, after being matched together, also using Pearson’s r to identify the best matches.
Statistical analysis
We used generalized linear mixed models (GLMM)45 to compare the velocity and the vertical amplitude of the CoMPACS, the length and duration of the gait cycles, the recovery rate, the external, internal, and total mechanical power produced by the PACS, the TotMCF and ∆MCF, the number of synergies, as well as the Pearson’s r-values across conditions.
The experimental condition was entered as a fixed factor in the model, and individuals as a random variable. We used a Gaussian GLMM for all variables, except for the comparison of the number of synergies across conditions, for which a Poisson GLMM was used. For Gaussian GLMMs we systematically inspected the normality of the model residuals with Q-Q plots. We used the functions lmer() and glmer() of the R package lme4 46 to run the Gaussian and the Poisson mixed models, respectively. The effect of experimental conditions was tested by comparing the deviance of the model with and without the fixed factor with a χ2 test. Multiple comparisons across experimental conditions were performed with the function glht() of the multcomp R package47 using the default Tukey test as post-hoc test. Pearson’s r were Fisher Z-transformed before running the analyses. The significance threshold was set to 0.05. All data in the text are given as mean ± SD. Since our sample size was low, which could lead to inflate the Type II error (not rejecting H0 when H0 is false), we followed the recommendations of Nakagawa & Foster (2004)49 and provide in the Supplemental Table S1 the value of Cohen d standardized effect size50, along with its 95% confidence interval51, for each studied parameter and each paired comparison between conditions. A confidence interval that largely extends on both sides of zero indicates an absence of effect that would probably not change with increasing the sample size.
Ethics statement
All methods used in this study were carried out in accordance with relevant guidelines and regulations. More

1.
Tilman, D., Balzer, C., Hill, J. & Befort, B. L. Global food demand and the sustainable intensification of agriculture. Proc. Natl Acad. Sci. USA 108, 20260–20264 (2011).
ADS  CAS  PubMed  Article  Google Scholar 
2.
Alexandratos, N. & Bruinsma, J. World Agriculture Towards 2030/2050: The 2012 Revision (FAO, ESA Working paper, Rome, 2012).

3.
Carlson, K. M. et al. Greenhouse gas emissions intensity of global croplands. Nat. Clim. Change 7, 63–68 (2017).
ADS  CAS  Article  Google Scholar 

4.
Houghton, R. A. et al. Carbon emissions from land use and land-cover change. Biogeosciences 9, 5125–5142 (2012).
ADS  CAS  Article  Google Scholar 

5.
Cassman, K. G., Dobermann, A., Walters, D. T. & Yang, H. Meeting cereal demand while protecting natural resources and improving environmental quality. Annu. Rev. Environ. Resour. 28, 315–358 (2003).
Article  Google Scholar 

6.
Laurance, W. F., Sayer, J. & Cassman, K. G. Agricultural expansion and its impacts on tropical nature. Trends Ecol. Evol. 29, 107–116 (2014).
PubMed  Article  Google Scholar 

7.
Lobell, D. B., Schlenker, W. & Costa-Roberts, J. Climate trends and global crop production since 1980. Science 333, 616–620 (2011).
ADS  CAS  Article  Google Scholar 

8.
Asseng, S. et al. Rising temperatures reduce global wheat production. Nat. Clim. Change 5, 143–147 (2015).
ADS  Article  Google Scholar 

9.
Rosenzweig, C. et al. Assessing agricultural risks of climate change in the 21st century in a global gridded crop model intercomparison. Proc. Natl Acad. Sci. USA 111, 3268–3273 (2014).
ADS  CAS  Article  Google Scholar 

10.
Mueller, N. D. et al. Closing yield gaps through nutrient and water management. Nature 490, 254–257 (2012).
ADS  CAS  PubMed  Article  Google Scholar 

11.
Jägermeyr, J., Pastor, A., Biemans, H. & Gerten, D. Reconciling irrigated food production with environmental flows for Sustainable Development Goals implementation. Nat. Commun. 8, 15900 (2017).
ADS  PubMed  PubMed Central  Article  CAS  Google Scholar 

12.
Schauberger, B. et al. Consistent negative response of US crops to high temperatures in observations and crop models. Nat. Commun. 8, 13931 (2017).
ADS  CAS  PubMed  PubMed Central  Article  Google Scholar 

13.
Tack, J., Barkley, A. & Hendricks, N. Irrigation offsets wheat yield reductions from warming temperatures. Environ. Res. Lett. 12, 114027 (2017).
ADS  Article  Google Scholar 

14.
Troy, T. J., Kipgen, C. & Pal, I. The impact of climate extremes and irrigation on US crop yields. Environ. Res. Lett. 10, 054013 (2015).
ADS  Article  Google Scholar 

15.
Elliott, J. et al. Constraints and potentials of future irrigation water availability on agricultural production under climate change. Proc. Natl Acad. Sci. USA 111, 3239–3244 (2014).
ADS  CAS  PubMed  Article  PubMed Central  Google Scholar 

16.
Jägermeyr, J. et al. Integrated crop water management might sustainably halve the global food gap. Environ. Res. Lett. 11, 025002 (2016).
ADS  Article  Google Scholar 

17.
Rosegrant, M. W., Ringler, C. & Zhu, T. J. Water for agriculture: maintaining food security under growing scarcity. Annu. Rev. Environ. Resour. 34, 205–222 (2009).
Article  Google Scholar 

18.
Siebert, S. & Döll, P. Quantifying blue and green virtual water contents in global crop production as well as potential production losses without irrigation. J. Hydrol. 384, 198–217 (2010).
ADS  Article  Google Scholar 

19.
Li, X. & Troy, T. J. Changes in rainfed and irrigated crop yield response to climate in the western US. Environ. Res. Lett. 13, 064031 (2018).
ADS  Article  Google Scholar 

20.
Neverre, N., Dumas, P., Nassopoulos, H. Large-scale water scarcity assessment under global changes: insights from a hydroeconomic framework. Hydrol. Earth Syst. Sci. Discuss. Preprint at https://doi.org/10.5194/hess-2015-502 (2016).

21.
Lesk, C., Rowhani, P. & Ramankutty, N. Influence of extreme weather disasters on global crop production. Nature 529, 84–87 (2016).
ADS  CAS  PubMed  Article  Google Scholar 

22.
Lobell, D. B., Sibley, A. & Ortiz-Monasterio, J. I. Extreme heat effects on wheat senescence in India. Nat. Clim. Change 2, 186–189 (2012).
ADS  Article  Google Scholar 

23.
Gourdji, S. M., Sibley, A. M. & Lobell, D. B. Global crop exposure to critical high temperatures in the reproductive period: historical trends and future projections. Environ. Res. Lett. 8, 024041 (2013).
ADS  Article  Google Scholar 

24.
Espe, M. B. et al. Point stresses during reproductive stage rather than warming seasonal temperature determine yield in temperate rice. Glob. Change Biol. 23, 4386–4395 (2017).
ADS  Article  Google Scholar 

25.
Martre, P. et al. Multimodel ensembles of wheat growth: many models are better than one. Glob. Change Biol. 21, 911–925 (2015).
ADS  Article  Google Scholar 

26.
Li, T. et al. Uncertainties in predicting rice yield by current crop models under a wide range of climatic conditions. Glob. Change Biol. 21, 1328–1341 (2015).
ADS  CAS  Article  Google Scholar 

27.
Müller, C. et al. Global gridded crop model evaluation: benchmarking, skills, deficiencies and implications. Geosci. Model Dev. 10, 1403–1422 (2017).
ADS  Article  Google Scholar 

28.
Wang, X. et al. Emergent constraint on crop yield response to warmer temperature from field experiments. Nat. Sustain. https://doi.org/10.1038/s41893-020-0569-7 (2020).

29.
Folberth, C. et al. Uncertainty in soil data can outweigh climate impact signals in global crop yield simulations. Nat. Commun. 7 https://doi.org/10.1038/ncomms11872 (2016).

30.
Elliott, J. et al. The Global Gridded Crop Model Intercomparison: data and modeling protocols for Phase 1 (v1.0). Geosci. Model Dev. 8, 261–277 (2015).
ADS  Article  Google Scholar 

31.
Roberts, M. J. et al. Comparing and combining process-based crop models and statistical models with some implications for climate change. Environ. Res. Lett. 12, 095010 (2017).
ADS  Article  Google Scholar 

32.
Raftery, A. E., Gneiting, T., Balabdaoui, F. & Polakowski, M. Using Bayesian model averaging to calibrate forecast ensembles. Mon. Weather Rev. 133, 1155–1174 (2005).
ADS  Article  Google Scholar 

33.
Monfreda, C., Ramankutty, N. & Foley, J. A. Farming the planet:2. Geographic distribution of crop areas, yields, physiological types, and net primary production in the year 2000. Glob. Biogeochem. Cycles 22, GB1022 (2008).
ADS  Article  CAS  Google Scholar 

34.
Fekete, B. M., Vorosmarty, C. J., & Grabs, W. Global, composite runoff fields based on observed river discharge and simulated water balances. Report No. 22. (World Meteorological Organization–Global Runoff Data Center, Koblenz, Germany, 1999).

35.
Wang, X. X. et al. Taking account of governance: Implications for land-use dynamics, food prices, and trade patterns. Ecol. Econ. 122, 12–24 (2016).
Article  Google Scholar 

36.
Neumann, K. et al. Exploring global irrigation patterns: a multilevel modelling approach. Agric. Syst. 104, 703–713 (2011).
Article  Google Scholar 

37.
Fekete, B. M. & Vörösmarty, C. J. The current status of global river discharge monitoring and potential new technologies complementing traditional discharge measurements. IAHS Publ. 309, 129–136 (2007).
Google Scholar 

38.
Hanasaki, N., Yoshikawa, S., Pokhrel, Y. & Kanae, S. A global hydrological simulation to specify the sources of water used by humans. Hydrol. Earth Syst. Sci. 22, 789–817 (2018).
ADS  Article  Google Scholar 

39.
Siebert, S., Henrich, V., Frenken, K., & Burke, J. Update of the digital global map of irrigation areas to version 5 (Rheinische Friedrich-Wilhelms-Universität, Bonn, Germany, Food and Agriculture Organization of the United Nations, Rome, Italy, 2013).

40.
Rosa, L. et al. Potential for sustainable irrigation expansion in a 3 °C warmer climate. Proc. Natl Acad. Sci. USA 117, 29526–29534 (2020).
ADS  CAS  PubMed  Article  Google Scholar 

41.
Rosa L., Chiarelli D. D., Rulli M. C., Dell’Angelo J. & D’Odorico P. Global agricultural economic water scarcity. Sci. Adv. 6, eaaz6031 (2020).

42.
Rosegrant, M. W. & Cai, X. Global water demand and supply projections. Water Int. 27, 170–182 (2002).
Article  Google Scholar 

43.
Berkoff, J. China: the South–North Water Transfer Project—is it justified? Water Policy 5, 1–28 (2003).
Article  Google Scholar 

44.
Schmitz, C. et al. Blue water scarcity and the economic impacts of future agricultural trade and demand. Water Resour. Res. 49, 3601–3617 (2013).
ADS  Article  Google Scholar 

45.
Biewald, A., Rolinski, S., Lotze-Campen, H., Schmitz, C. & Dietrich, J. P. Valuing the impact of trade on local blue water. Ecol. Econ. 101, 43–53 (2014).
Article  Google Scholar 

46.
Portmann, F. T., Siebert, S. & Döll, P. MIRCA2000—Global monthly irrigated and rainfed crop areas around the year 2000: a new high-resolution data set for agricultural and hydrological modeling. Glob. Biogeochem. Cycles 24, GB1011 (2010).
ADS  Article  CAS  Google Scholar 

47.
Muller, C. et al. The Global Gridded Crop Model Intercomparison phase 1 simulation dataset. Sci. Data 6, 50 (2019).
PubMed  PubMed Central  Article  Google Scholar 

48.
Kobayashi, K. & Salam, M. U. Comparing simulated and measured values using mean squared deviation and its components. Agron. J. 92, 345–352 (2000).
Article  Google Scholar 

49.
Deutscher Wetterdienst Frankfurt M. Mitteilungen des Deutschen Wetterdienstes: Windschutzanlagen auf der hohen Rhön (Dt. Wetterdienst P.11, 1954).

50.
Gornott, C. & Wechsung, F. Statistical regression models for assessing climate impacts on crop yields: a validation study for winter wheat and silage maize in Germany. Agric. Meteorol. 217, 89–100 (2016).
Article  Google Scholar 

51.
Harris, I., Jones, P. D., Osborn, T. J. & Lister, D. H. Updated high-resolution grids of monthly climatic observations – the CRU TS3.10 Dataset. Int J. Climatol. 34, 623–642 (2014).
Article  Google Scholar 

52.
Jägermeyr, J. et al. Water savings potentials of irrigation systems: global simulation of processes and linkages. Hydrol. Earth Syst. Sci. 19, 3073–3091 (2015).
ADS  Article  Google Scholar 

53.
Pastor, A. V., Ludwig, F., Biemans, H., Hoff, H. & Kabat, P. Accounting for environmental flow requirements in global water assessments. Hydrol. Earth Syst. Sci. 18, 5041–5059 (2014).
ADS  Article  Google Scholar 

54.
Poff, N. L. et al. The natural flow regime: a paradigm for river conservation and restoration. BioScience 47, 769–784 (1997).
Article  Google Scholar  More