Philippa Kaur

Study site and future climate treatmentOur study site is located at the Bad Lauchstädt Field Research Station, Bad Lauchstädt, Germany (51° 22060 N, 11° 50060 E), which belongs to the Helmholtz Centre for Environmental Research–UFZ. Long-term mean annual precipitation in the area is 489 mm and the mean annual temperature is 8.9 °C (ref. 32). During 2018 and 2019, Europe experienced a record-setting drought that was especially severe in 2018 (refs. 33,34); the mean annual precipitation at our study site in 2018 and 2019 was 254 mm and 353 mm, respectively, whereas 2017 was a more normal year, with a mean annual precipitation of 403 mm. Mean annual temperatures were above average: 2017, 10.5 °C; 2018, 10.8 °C; 2019, 11.2 °C (data from the weather station at the Bad Lauchstädt field station). The soils in the study area are fertile Haplic Chernozem type32,35.Our eDiValo experiment was conducted in the GCEF, which was designed to investigate climate change effects under different land-use scenarios32. We used 10 ‘extensively’ used pastures of the GCEF in our experiment; that is, 384-m2 (16 × 24 m) areas of grassland (hereafter called ‘pastures’) that were grazed by a flock of 20 sheep 2–3 times each year. Grazing was implemented as short-time high-intensity grazing events, each lasting 24 h (ref. 32). This type of high-intensity but short-term grazing is considered better in maintaining species richness as it gives plants more time to recover between grazing events36. It is also a recommended management type for nature conservation areas in Germany37. Vegetation in the pastures was species-rich grassland vegetation that is typical of drier regions of central Germany32,38. The whole GCEF was fenced to exclude native large mammalian herbivores (for example, deer); however, European hare (Lepus europaeus), wood mice (Apodemus sylvaticus) and voles (Microtus arvalis) are common at the site.Our experimental design was originally intended to test the dependence of light competition on nutrient and herbivory under current and future climatic scenarios. Although we included both climate treatments in our data, climate was never significant for richness and Shannon diversity, either alone or in interaction with other factors, and our focus was therefore on the other treatments. Five of the above random pastures received future climatic treatment which was based on different dynamic regional climate models for Germany, all predicting an increased mean temperature by approximately 2 °C year-round, strongly decreased summer precipitation and slightly increased spring and autumn precipitation (https://www.regionaler-klimaatlas.de/) (ref. 32). Passive night-time (after sunset and before sunrise) warming through the use of roller blinds attached to the GCEF roof and eastern and western wall structures was used to increase the air temperature. In each spring (1 March–31 May) and autumn (1 September–30 November), future climate plots received 110% of the ambient rainfall and in the summer (1 June–31 August), they received 80% of the ambient rainfall. The precipitation treatment was adjusted weekly and compensated for a possible night-time reduction in rainfall due to temperature treatment. A detailed description of the future climate treatment is provided in a previous report32.Fertilization, herbivore exclusion and light additionWe first tested whether adding light can offset the negative effect of fertilization on plant diversity. In May 2017, we established a full-factorial experiment of fertilization and light addition. Within each 10 pastures (5 in ambient climatic conditions, 5 in future climatic conditions), we established 4 plots of 1.4 × 1.4 m, separated by a 1-m buffer zone (hereafter called ‘blocks’), in total 40 plots and 10 blocks. At the time the experiment was established, vegetation in the whole experimental area (that is, in a block of 4 plots and the surrounding 1-m area) was trimmed to a height of 5 cm to make conditions uniform and the whole area was temporarily fenced to let the experiment establish and fertilization effects develop. The temporary fence was removed in August when the herbivore exclusion treatment was started. Therefore, there was no grazing by sheep in the experimental plots in the summer of 2017. Two randomly chosen plots received fertilizer treatment and two were controls. For the former (fertilizer-treatment plots), slow-release granular NPK fertilizer (a mixture of Haifa Multicote 2 M 40-0-0 40% N; Triple Super Phosphate (TSP) 45% P205; and potassium sulfate fertilizer 50% K2O, 45% SO3) was added twice per growing season, in a total of 10 g N, 10 g P and 10 g K per m² (see ref. 3 for a similar protocol that is used in grasslands worldwide). In 2017, the first fertilization was done at the beginning of June right after establishing the experiment and the second fertilization was done at the beginning of July. In the subsequent years, the first fertilization was done at the beginning of the growing season (late March–April) and the second fertilization was done in June. In 2019, two previously unfertilized plots were accidentally fertilized and were thereafter treated as fertilized plots. To manipulate light, 1.4 × 1.4-m plots were further divided into two subplots, 0.7 m × 1.4 m each, and one of these was randomly assigned to the light-addition treatment, resulting in 80 subplots (Fig. 1). We installed two 120-cm-long and 3.5-cm-wide recently developed LED lamps (C65, Valoya) parallel to each other and at a 28-cm distance from each other to each light-addition subplot. To increase light for the small understory plants that are the most likely to suffer from competition for light, we installed the lamps 10 cm above the smallest plants. The lamps were gradually uplifted over the course of the growing season to follow the growth of the smallest plants. As our light-addition treatment was intended to mimic natural sunlight (that is, making a gap in a dense vegetation and allowing the sunshine in), we chose the spectrum of the lamps to include all wavelengths of sunlight, including small amounts of ultraviolet and infrared. Each lamp added roughly 350–400 µmol and did not alter the air or aboveground soil surface temperature (Fig. 1b), which is an improvement on previous studies12. Each year, we added light during the active growing season: the lamps were switched on early in the spring (March–April), when temperatures were clearly above zero, and switched off and removed when temperatures dropped close to zero in November–December and aboveground plant parts had died and formed litter. Each day, the lamps were set to switch on two hours after sunrise, and to switch off two hours before sunset, and when the temperature exceeded 28 °C to prevent overheating. We did not install unpowered lamps to unlighted plots because our modern, narrow LED lamps caused minimal disturbance (see below) and no heating (Fig. 1b), and because unpowered lamps would have added an artefact in that they create shade that does not occur when the lamps are on in lighted plots.At the end of August 2017, after running the fertilization–light-addition experiment for one growing season, we expanded the experiment by implementing the herbivore exclusion treatment in a full-factorial combination with the other treatments. Two of the previously established 1.4 m × 1.4-m plots, one with and one without the fertilization treatment, were randomly allotted to the herbivore (sheep) exclusion treatment and fenced with rectangular metal fences of 1.8 m × 1.8 m, 82 cm height and 10 cm mesh size. At the same time, the temporary fence established in May 2017 was removed from around the whole experimental area, allowing the grazing of sheep in unfenced plots. The fences did not exclude mice, voles and hares. For the time of each grazing event, lamps in grazed subplots were removed and switched off in the ungrazed subplots. Uplifting the lamps from grazed plots did not cause disturbance because vegetation in grazed plots was always short and did not reach above the lamps. Inside exclosures, lamps were always kept in place during the growing season, and plants could freely grow around and above them.Plant community and trait samplingIn July 2017, we established 50 cm × 50-cm permanent quadrats in every subplot for plant community sampling. We visually estimated the per cent areal cover for all species occurring in the quadrats, and litter cover, from the beginning of June to mid-June 2019, when the vegetation was at its peak biomass. The 2017 sampling happened later, in mid-July, because vegetation in all plots and surrounding areas was trimmed to a height of 5 cm at the time of the establishment of the experiment at the end of May, and it took later for vegetation to reach its peak biomass. In 2018, the effects of drought were devastating, and most plants had senesced or died before the planned sampling date; we therefore omitted the year 2018. At the beginning of each growing season—that is, when the lamps were installed and switched on—there was very little live biomass in the plots, and the maximum height of existing plants was approximately 5 cm (in all plots). During the peak biomass the maximum plant height was up to approximately 1 m; however, it varied greatly between the treatments and was especially low in grazed plots. All vegetation surveys were done by the same trained and experienced person with a minimum estimate threshold of 0.1%. We used plant cover data to calculate species richness and Shannon diversity.In May–June 2020, we measured plant height (centimetres), SLA (leaf area in square millimetres per milligram of dry mass), foliar C:N (based on the per cent C and N in plant leaves) and LWC (leaf water content as 1,000 − LDMC (the ratio of leaf dry mass to saturated fresh mass), expressed as milligrams per gram39) for most species occurring in the experimental plots, and complemented the trait data from the TRY Plant Trait Database40,41,42 (v.5.0; https://www.try-db.org/TryWeb/Home.php) and for one species one trait value from another source9. The trait data were collected from seven to ten individuals per species from the study site or close areas; the collection and handling followed standard protocols39. We chose these traits because they are widely documented to be associated with responsiveness to soil nutrients, herbivory and light9,26,27,43,44,45,46. We used all traits as, although they partially reflect similar ecological adaptations (for example, leaf economics spectrum43), they could also potentially reflect independent and distinctive processes, and differently mediate the responses of species to our treatments. For example, SLA and LWC in our dataset correlated weakly (r2 = 0.16), but were to a greater extent uncorrelated (Extended Data Table 6), and could function differently, for example, in light capture and drought tolerance26,39. In 2017, our trait data covered on average 97.7–98.6% of the total cover in the plots, the value slightly differing depending on the trait as we did not have all traits for all species. Our own trait collections covered on average 96.6–97.6% and TRY data covered on average 0.9–2% of the total cover. In 2019, the whole trait data covered on average 99.5% of the total cover in the plots, again slightly depending on the trait. Our own trait collections covered on average 94.2–96.5% and TRY data covered on average 2.7–5.3% of the total cover.Abiotic environmental measurementsWe measured several soil and other environmental properties from the experimental plots. Light availability (photosynthetically active radiation; PAR) in unlighted and lighted (under lamps) subplots was measured using LI-190R and LI-250A meters (LI-COR), approximately 7–10 cm under the lamps and 15–20 cm above ground level. We measured light availability from the same distance to the ground in unlighted plots. Measurements of light availability were done in mid-July 2020 on three consecutive cloudless days around noon. Note that in grazed plots, light levels between lighted and unlighted plots are more similar than inside exclosures (Fig. 1), because herbivores keep the vegetation short, and natural sunlight can therefore reach under the lamps where the light measurements were taken. Air temperature and humidity were recorded from unlighted and lighted (under lamps) subplots using loggers (HOBO MX2301A, Onset Computer Cooperation) that were installed approximately 7 cm under the lamps and to the same height from the ground in unlighted plots, and were replicated under different combinations of fertilization, herbivore exclusion and light addition in ambient climatic conditions three times (n = 3). The logger data were collected in May 2019 before the effects of drought were visible.Statistical analysisWe analysed our data in two steps. First, to test whether competition for light mediates the effect of fertilization on diversity, we analysed the effects of fertilization and light and their interaction on species richness and Shannon diversity using data from 2017, when the herbivore exclusion treatment had not yet been implemented. We also analysed the effects of treatment on total vegetation cover and litter cover. We fit LME models in which diversity (species richness and Shannon diversity), total cover and litter cover, each in their own model, were explained by fertilization, light addition and their interaction (fixed variables). All treatments were categorical variables with two levels (treated and untreated). In each model, subplot was nested within plot, which was nested within block (nested random variable). We simplified the models using the anova() function for model comparison in the nlme and lme4 packages in R (ref. 47) (on the basis of log likelihood ratio tests; P ≥ 0.05; Extended Data Table 2). This was done to uncover the significance of the main effects and interaction terms, to avoid overparametrization47,48 and to provide model-derived parameter estimates for the figures (Extended Data Table 5). However, we also provide full model results that are qualitatively similar to the results of simplified models (Extended Data Tables 3 and 4); therefore, model choice did not affect our conclusions. Climate treatment was included in all original models but was never significant for richness and diversity, and was not considered further. Total cover and litter results for 2017 are reported in Extended Data Figs. 1a,b and 3a). As there was heterogeneity in the variance structure between treatments, we used the varIdent() function in the nlme package in R to allow each treatment combination to have a different variance. Model fit was inspected using model diagnostic plots in the package nlme. In the full design with climate included, the number of replicates per treatment combination was ten.Second, to include herbivore exclusion to the experimental design and to test whether competition for light mediates the effect of herbivore exclusion on diversity, and whether competition for light, herbivory and fertilization interact, we analysed the effects of herbivore exclusion, fertilization, light and their interactions on species richness and Shannon diversity using data from 2019. All treatments were categorical variables with two levels (treated and untreated). We also analysed the effects of treatment on total vegetation cover and litter cover. We fit similar models to those described above, except that herbivore exclusion was an additional fixed factor in the models. We simplified the models, used the varIdent() function to account for heteroscedasticity and checked the model fit using model diagnostic plots, as above. Climate treatment was included in all original models but was significant for litter cover only, and was not considered further. In the full design with climate included, the number of replicates per treatment combination was five.To further assess which plant traits increased the probability of species benefiting from the addition of light, we first created a binary response variable: those species that increased from unlighted to lighted plots (that is, had a higher value in a lighted than an unlighted plot) were given a value of 1 and those that did not were given a value of 0. This response variable takes into account rare species that emerged or persisted in the lighted plots but were absent in the unlighted plots (that is, species gains and losses) and changes in small, subordinate species (those that are likely to benefit from light addition) with small but consistently trait-dependent changes in response to light. It is also in line with our species richness analyses, as species gains and losses ultimately determine richness responses. We did not use different indexes (for example, lnRR or RII) because these could not handle multiple zero values and species losses or gains (that is, species having zero cover in either unlighted or lighted subplots). Second, we fit GLME models with a binomial error structure (family = “binomial”, link = “logit”) in which a probability of a species increasing from unlighted to lighted plots was explained by categorical experimental treatments (fertilization, herbivore exclusion and their interactions), traits (SLA, height, LWC, foliar C:N), and interactions between the treatments and traits. Each trait was analysed in its own model as some of the traits were correlated (Extended Data Table 6), and to avoid overly complex models and overparametrization47,48. We included all species for which we had traits in the models. As we calculated the increase in cover from unlighted to lighted plots, our smallest experimental unit in trait analyses was a plot (not a subplot, unlike in other analyses). As there were several species in the same plots, we nested species within plots, and plots within blocks. We similarly simplified the models to include only significant variables (on the basis of χ2 tests; P ≥ 0.05). We did not include a crossed random effect for species in the models because the full models with a more complex random structure did not converge; however, when we refitted the simplified models with a crossed random effect for species, we found that the models converged (with scaled data) and that the significance of the effects remained qualitatively the same. Climate was included in all original models but was never significant. In addition, C:N and height did not predict the responsiveness of species to light in either year (P ≥ 0.13 for both); results are therefore not shown. In the full design with climate included, the number of replicates per treatment combination was five; however, the number of observations was greater (see Fig. 4 and Extended Data Fig. 4). To make sure that our results for SLA and LWC were not influenced by whether they were analysed in separate models or in the same model, or by the order in which they were in the models, we also performed analyses in which both SLA and LWC were included (in both orders). Results remained qualitatively similar and are not discussed further.Furthermore, to check whether our trait results were driven primarily by species gains and losses or changes in abundance, we ran additional trait analyses for which we calculated the change in cover between lighted and unlighted subplots (cover in lighted subplot − cover in unlighted subplot), and analysed the ‘change’ with otherwise similar trait models to those described above, except that we used Gaussian error structure. With this index, which gives a disproportionate importance to the abundant species, we found that traits were poor predictors of changes in cover between lighted and unlighted plots (all interactions were non-significant, P  > 0.05, except for a marginally significant C:N × fertilization interaction in 2017 that was no longer visible in 2019; results not shown; codes and data available in the Dryad repository). We also analysed presence–absence-based species losses and gains. In these models, each species was given a value of 1 when it was present in the lighted subplot but absent from the unlighted subplot; otherwise, these models were similar to the binomial trait models described above. These models produced, to a large extent, similar results to our models using the probability of increase in response to light as a response variable (results not shown; codes and data available in the Dryad repository). These additional analyses and results support using the probability of increase in response to light as our response variable, rather than abundance-based metrics, as it includes both gains and losses and abundance aspects, and is therefore a general test that is well suited to assessing species gains and extinctions and changes in subordinate species.All statistical analyses were performed using R v. 4.0.0 (ref. 49). We used the nlme package (v.3.1.147) for LME models50, the lme4 package (v.1.1.23) for GLME models51, and the car package52 for P values (v.3.07).Reporting summaryFurther information on research design is available in the Nature Research Reporting Summary linked to this article. More

The link between SDGs and NRRPThe Italian National Recovery and Resilience Plan (NRRP) is part of the Next Generation EU (NGEU) program, the 750-billion-euro package, consisting of about half of grants, agreed by the European Union in response to the pandemic crisis. The main component of the NGEU program is the Recovery and Resilience Facility (RRF), which has a duration of six years, from 2021 to 2026, and a total size of €672.5 billion (€312.5 billion grants, the remaining €360 billion loans at subsidized rates).The Plan is developed around three strategic axes shared at European level: digitalization and innovation, ecological planning and social inclusion.The missions of the NRRP are as follows:

Mission 1: Digitalization, innovation, competitiveness, culture and tourism

Mission 2: Green revolution and ecological transition

Mission 3: Infrastructure for sustainable mobility

Mission 4: Education and research

Mission 5: Cohesion and inclusion

Mission 6: Health.

With the aim of encouraging the debate on the use of sustainability indicators for monitoring the progress of the PNRR, a mapping of the correspondences between the 17 Sustainable Development Goals and the 6 Missions provided for by the NRRP is proposed (Fig. 1). In this way it is possible to identify the SDGs indicators that can be useful tools for achieving the missions of the NRRP.Figure 1Relationships between SDGs indicators and NRRP missions.Full size imageOf particular interest for the purposes of our work is Mission 2 (Green Revolution and Ecological Transition) of NRRP. It provides for investments and reforms for the circular economy and to improve waste management, strengthen separate collection infrastructure and modernize or develop new waste treatment plants. Substantial tax incentives are provided to increase the energy efficiency of buildings, to achieve progressive decarbonization, to increase the use of renewable energy sources. In addition, the Mission devotes resources to enhancing the capacity of electricity grids, their reliability, security, and flexibility (Smart Grid) and water infrastructure. The Mission also includes the issues of territorial security, with prevention and restoration interventions in the face of significant hydrogeological risks, the protection of green areas and biodiversity, and those related to the elimination of water and soil pollution, and the availability of water resources.The main components of this mission are:

M2C1: Circular economy and sustainable agriculture

M2C2: Renewable energy, hydrogen, grid, and sustainable mobility

M2C3: Energy efficiency and upgrading of buildings

M2C4: Protection of land and water resources.

The analysis of Mission 2 (Green Revolution and Ecological Transition) finds ample space in the SDGs creating important interconnections between the different indicators present in the individual Goals and the objectives of the Mission itself.The SDGs indicators to support the NRRPThe SDGs indicators selected for the analysis of Mission 2 (Green Revolution and Ecological Transition) of the NRRP, are descripted in Table 1. We considered 13 indicators, selected from Goals 2, 6, 7, 11, 12 and 15 which may be of significant interest for the achievement of Mission 2. These indicators will then be attributed to the individual components of the mission.Table 1 Goal, indicators, measures e source of SDGs data.Full size tableThe indicators were chosen based on their relevance to the objectives of the mission and on the availability of data on a regional basis. For each main component we can use the following indicators:

M2C1: Circular economy and sustainable agriculture:

– Share of utilized agricultural area invested by organic crops

– Growth rate of organic crops

– Delivery of municipal waste to landfill.

– Separate waste collection

M2C2: Renewable energy, hydrogen, grid and sustainable mobility:

M2C3: Energy efficiency and upgrading of buildings

M2C4: Protection of land and water resources

– Irregularities in water distribution

– Sealing and soil consumption per capita

– Soil sealing from artificial cover

– Fragmentation of the natural and agricultural territory

– Incidence of urban green areas on the urbanized surface of cities.

The SDGs indicators at the level of territorial distribution in ItalyWe carry out a first analysis by territorial distribution for the different sets of main components of Mission 2.From a first analysis of the M2C1 indicators (Circular Economy and Sustainable Agriculture) it emerges that the share of agricultural area destined for organic crops is greater, especially in the Center and in the South of Italy. In 2019, the extent of organic farming in Italy reached 15.8% of the utilized agricultural area, almost double the EU average. However, the annual growth rate of the areas converted to organic farming or in the process of conversion (+ 1.8%) is the lowest since 2012 and is negative in the South, where for the second consecutive year there is a decrease (− 2.1% in the 2-year period 2017–2019). The dynamics of organic farming is an index of the spread of sustainable agricultural practices, which must be accompanied by measures that also consider the pressure on the environment generated by agriculture (Table 2).Table 2 M2C1 indicators—Circular economy and sustainable agriculture by territorial distribution (year 2019).Full size tableAlso, in the Central and Southern Italy area there is the greatest delivery of waste to landfills. Waste cycle management is crucial for living conditions and global health. The share of municipal waste landfilled is steadily decreasing at national level. In 2019, in fact, the part sent to landfill is equal to 20.9% of the total, down compared to the previous year (21.5%). The separate collection of municipal waste represents a further important step in view of the objective of reducing the amount of waste returned to the environment and, more specifically, of the delivery of waste to landfills. The 18.5 million tons of differentiated RU in 2019 represent 61.3% of national production, a share almost doubled compared to ten years ago and up from last year by 3.1 percentage points. Despite the evident progress, Italy is still marked by a considerable delay compared to the regulatory objectives, having not yet reached, in 2019, the target of 65% of separate collection planned for 2012. Critical issues are also observed in relation to the substantial territorial gaps, which disadvantage the Center and the South compared to the North, despite the distances have been reduced in recent years.
Regarding the M2C2 Mission (Renewable Energy, Hydrogen, Network and Sustainable Mobility), national and international energy policies have been committed for years to the enhancement of renewable energy sources, with the aim of decarbonizing the economy and guaranteeing the commitments made in the field of climate change. In 2019, one year after the expiry of the objectives of the European Union’s Climate-Energy Package, fourteen Member States, including Italy, exceeded the target assigned at national level. In Italy, the overall share of energy from renewable sources in gross final consumption (CFL) of energy, equal to 18.2% (Table 3), a percentage slightly lower than the average of the EU27 (19.7%), is placed for the sixth consecutive year above the 17% target set for our country. However, for Italy to achieve the ambitious programs defined by the 2020 National Integrated Energy and Climate Plan, which set a 30% target for renewables by 2030, a further boost to production from renewable sources is necessary. The resources introduced by the National Recovery and Resilience Plan (NRRP) to achieve the “green revolution and ecological transaction” include significant investments in the energy field, focusing, among other components, on a further strengthening of the Sources from Renewable energy (FER).Table 3 M2C2 indicators—Renewable energy, hydrogen, network and sustainable mobility by territorial distribution (year 2019).Full size tableThe M2C3 Mission (Energy Efficiency and Upgrading of Buildings) devotes resources to enhancing the capacity of electricity grids, their reliability, safety, and flexibility (Smart Grid). Consistent with the objectives of reducing energy consumption pursued by European policies, the Italian figure for 2019 confirms the process of reducing Italian energy intensity, which marks a further contraction of 1.3%, reaching an overall negative balance compared to the last decade of 11.8%, with an average annual rate of change of − 1.2% (Table 4). The reduction in energy intensity is largely attributable to the effect of the measures in favor of energy efficiency, which, between 2011 and 2019, resulted in energy savings of 12 Mtoe/year, equal to 77% of the 2020 target set by the National Action Plan for Energy Efficiency 2017. A further acceleration of energy efficiency is expected, in the coming years, because of the investment plan envisaged by the NRRP, also linked to the redevelopment of the public and private building stock. At the sectoral level, the reduction in energy intensity is driven by improvements in industry, which, despite the slight increase in the last year, in 2019, with 92 toes per million euros, shows a decrease compared to 2009 of 17%, with an average annual rate of change of − 1.8%.Table 4 M2C3 indicators—Energy efficiency and requalification of buildings by territorial distribution (year 2019).Full size tableThe M2C4 Mission (Protection of the territory and water resources) also includes the issues of territorial safety, with prevention and recovery interventions, the protection of green areas and those related to the elimination of water and soil pollution.Italy is among the European countries of the Mediterranean area that use groundwater, springs and wells the most; these represent the most important resource of fresh water for drinking water use on the Italian territory (84.8% of the total withdrawn). The efficiency of municipal drinking water distribution networks has been steadily deteriorating since 2008 for more than half of the regions. The share of families who complain of irregularities in the water supply service in their home is stable (equal to 8.6% in 2019) with more accentuated values in the Center and South of Italy (Table 5).Table 5 M2C4 indicators—Protection of land and water resources by territorial distribution (year 2019).Full size tableLand degradation, understood as loss of ecological functionality, is monitored through the dynamics of land consumption, which Italy has committed to zero by 2030 with the National Strategy for Sustainable Development (2017). The “consumed” soil is that occupied by urbanization and made impermeable by artificial roofing (soil sealing). Excessive fragmentation of open spaces, however, is also a factor of degradation, since the barriers made up of buildings and infrastructures interrupt the continuity of ecosystems, making even unoccupied but not large enough spaces ecologically inert and unproductive. Moreover, in a fragile territory such as Italy, land consumption is also a significant factor of hydrogeological risk and deterioration of the landscape. The index of sealing and land consumption per capita in 2019 increases for the fifth consecutive year, resulting in 357 m2 per inhabitant. The soil sealed by artificial covers is equal to 7.1% of the national territory (8.5% in the North, 6.7% in the Center, 5.9% in the South).According to Ispra estimates, 44.3% of Italy’s natural and agricultural land has a high or very high degree of fragmentation. A joint representation of the variations in fragmentation and soil sealing over the last two years summarizes recent trends in land consumption and their impact on the environment and landscape.A further objective for 2030 is to provide universal access to safe, inclusive, and accessible public green spaces, for women and children, the elderly, and people with disabilities. In 2019 the incidence of urban green areas on the urbanized surface of cities is equal to 8.5% in Italy with slightly higher values in the North and less elevated in the South. More

Birds and housing62 males and 8 females of house sparrows were caught with mist nets in September and October 2019 in several sites in Kraków, Poland. Before releasing them to the outdoor aviary located on the campus of the Jagiellonian University, Kraków, Poland, each bird was weighed and banded with a metal band. The aviary measured 3.5 m in width, 10.0 m in length, 2.5 m in height, and was outfitted with trees, bushes, perches, wooden shelters, a water source, and food dishes. Initially, birds were maintained with water and a mixture of seeds: wheat, barley, millet, and sunflower seeds, provided ad libitum. Additionally, they had access to sand with shells and sepia.Experimental designAfter a few weeks of acclimation to captivity, the aviary was divided into two separate parts (3.5 × 5 m): aviary no. 1 (A1) and aviary no. 2 (A2). At the same time male individuals were assigned to two crossed experimental treatments, ensuring that in each aviary birds originated from all sampled populations. The experiment comprised of two different treatments conducted simultaneously—one designed to simulate a deficiency in an environmental factor influencing colouration (the quality of available food), the other—to introduce physiological stress and facilitate trade-offs in the allocation of resources limited by the first treatment (an immune response induced by a bacteria-derived compound, S1).The dietary manipulation was achieved by feeding one group of birds with a low-quality protein food (diet reduced in exogenous amino acids, namely phenylalanine and tyrosine content, which are precursors essential for melanin synthesis; PT-reduced diet), and the other one with a wholesome diet (control diet). At the same time, two levels of immune challenge were achieved within each dietary group, by injecting half of the birds with either lipopolysaccharide (LPS) from the cell wall of Escherichia coli, or a 0.9% saline vehicle (as a control). Four females were placed in each group of males to alleviate interspecific conflicts occurring in all-male sparrow flocks, but they did not take part in the experiments. After three weeks of experiment, birds housed in A1 were moved to A2, whereas birds from A2 were moved to A1.Immune challengeBefore receiving injections, birds were first weighed and then transferred from the outdoor aviary to the laboratory. 31 house sparrows (from both dietary groups) were injected intraperitoneally with 0.026 mg LPS (serotype O55:B5, Sigma-Aldrich) diluted in 0.1 mL of 0.9% saline vehicle, so that each bird received a dose of ca. 1 mg/kg body mass, which had previously been shown to induce sickness behaviour in another passerine, the white-crowned sparrow, Zonotrichia leucophrys55. 31 control males were injected with the same volume (0.1 mL) of 0.9% saline vehicle. All individuals were injected twice throughout the experiment with an interval of three weeks between the injections. Birds were always injected at the same time in the morning and early afternoon (between 9:00 am and 12:30 pm).Diet manipulationDuring the six weeks of the experiment (S1), birds received synthetic diet ad libitum, which constituted of a mixture of protein (WPC80, free amino acids and whey protein isolate BiPRO GMP 9000 (Agropur Inc., Appleton, USA)), fats, carbohydrates, and fiber30. The ingredients were thoroughly mixed to produce small pellets (6 mm in diameter) that the sparrows consumed readily. The experimental diet had phenylalanine and tyrosine at 42% (N = 32) of their level in the control diet (N = 30)30. The food pellets were prepared by ZooLab (zoolab.pl/en/home, Sędziszów, Poland). Each bird was weighed before and after the experiment to monitor potential effects of diet on body mass of each animal. Following the experiment, during the next three consecutive days, the amount of food consumed by passerines within every 24 h (starting from 10 am each day to 10 am next day) was noted for both compartments of the aviary. Because of different numbers of individuals per aviary, an overall weight of food consumed in A1 and A2 was calculated per individual, respectively.Feathers samplingMoult of the black bib feathers was stimulated at the end of the moulting period occurring in natural conditions in early November. At day 1 of the dietary/immunological experiment (S1) a small area of the bib (around 25 mm2) was plucked from each male sparrow held in A1. At day 2 the same procedure was performed on individuals from A2. The time difference is orders of magnitude smaller than the timescale of feather growth and hence it would not affect the results in any way.Because the feather growth rate may differ during melanogenesis, with consequences for final colouration (if feathers grow at a faster rate, pigments may be deposited over a larger surface and therefore result in less intense colouration56, we measured the rate of feather development during the course of the experiment. After three weeks of the experiment, three feathers from the upper, central, and lower region of the previously plucked bib were plucked once again. The mass of the collected feathers was determined to the nearest 0.01 mg (XP26 Micro Balance, Mettler-Toledo, Greinfensee, Switzerland). The experiment was completed after six weeks after fully regrown and developed feathers from the bib and PC2 were sampled the second time (S1). Three feathers from the central part of previously plucked bib region were collected to perform transmission electron microscopy (TEM) imaging, whereas the feathers obtained from the rest of the regrown bib area were subjected to electron paramagnetic resonance (EPR) spectroscopy and feather microstructure analyses (greater spatial density of melanized barbs or barbules may affect colouration17.Feathers measurementsReflectance measurementsAn USB4000 spectrophotometer (range 300–700 nm) with the PX-2 Pulsed Xenon Lamp (Ocean Optics, Dunedin, FL, USA) and a bifurcated probe with 7 × 400 μm optical fibres, equipped with a permanently attached 3 mm long black collar, was used to quantify the brightness of the bib feathers collected at the end of the experiment. The measurements were taken with 90 ms integration time and the probe held at 90° to a feather’s surface. Calibration measurements of a Spectralon white standard (Ocean Optics. Largo, FL, USA) were taken every 15 min during measurements. The order in which the samples were measured was randomized in terms of belonging to the experimental group. From each sample (N = 62), seven feathers were chosen and stacked in one pile on a piece of black paper. Ten reflectance measurements were taken on each pile, avoiding distal, brighter parts of the feathers. The obtained spectra were averaged and smoothed in the package ‘pavo’57. Brightness was calculated as a sum of the reflectance values over all wavelengths of a spectrum, and its lower values were interpreted as those indicative of a more melanin-rich feathers (i.e., absorbing more light).Feather developmentEach feather (3 per individual; N = 62 individuals) was laid on a white card and covered by a microscope slide to flatten the naturally curved feathers. Digital photographs were taken using camera (Canon EOS 7D) and imported to ImageJ v1.52a Software (National Institutes of Health, USA). The lengths of fully developed and undeveloped (still in sheath) parts of each feather were measured. To estimate the degree of a feather’s development, the length of the developed part of the vane was divided by its total length (quill with rachis plus the developed vane, Fig. 4A).Figure 4House sparrow feathers sampled from bib after three weeks of the experiment. Feathers during development (A), a TEM cross-sections of feather sampled from bib after the experiment (B).Full size imageFeather densityBarb density measurements were performed on the sampled regrown black bib feathers (N = 2–3 for each individual; N = 62 individuals), but because of their sparser structure we calculated the number of non-down (i.e., rigid) barbs on both sides of the vane, and divided this number by two (to obtain an average single-sided number of barbs) and then by the length of the rachis.Melanosome density (TEM)Feathers sampled from the bib of male sparrows (N = 62) were fixed for transmission electron microscopy (TEM) analysis in a mixture of 0.25 M sodium hydroxide and 0.1% Tween for 20 to 30 min on a bench-top shaker. Next, the feathers were treated with formic acid and ethanol in the ratio of 2:3 for 2.5 h and dehydrated twice for 20 min in 100% ethanol. Samples were embedded in a mixture of the PolyBed 812 resin (20 ml), DDSA (9 ml), NMA (12 ml) and DMP-30 (0.82 ml). Resin infiltration was gradual from 15% resin content in ethanol through 50%, 70% to 100% without alcohol. Each step lasted for 24 h. Then, the feathers were placed in silicone embedding moulds (Agar Scientific) and transferred to an oven. The polymerization proceeded at the temperature of 60 °C for 16 h. The epoxy resin blocks were then trimmed to get rid of excess resin. The surface of each block was prepared by its trimming, starting from the end of the feather, to approximately 5 mm using a glass knife. Next, ultrathin sections (70 nm) were cut with a diamond knife (DIATOME A. G., Berno, Switzerland) on a microtome (UC7, Leica, Wetzlar, Germany) and collected on single slot grids coated with a formvar film. The sections were then contrasted in uranyl acetate and lead citrate for 3 min. They were viewed and photographed with a transmission electron microscope (TEM) JEOL 2100HT (Jeol Ltd, Tokyo, Japan) for the purpose of investigating the number and density of the embedded pigment granules. For each individual three photographs of the cross-sections from a similar feather region were selected. Melanosome density was measured as the number of melanin granules observed in the barb cross-section divided by its area. Images were analysed using Adobe Photoshop (cross-sections area) and ImageJ (number of melanosomes, Fig. 4B).Melanin content: electron paramagnetic resonance (EPR) spectroscopyQuality and quantity of melanin pigments58 in individual feather samples obtained from the bib of house sparrows (N = 57) were characterized using a Varian E3 spectrometer (Varian, Sunnyvale, LA, USA) equipped with a rectangular resonance (TE 102) cavity. Five milligrams of feathers per individual were placed inside the Wilmad finger quartz dewar WG-816-Q (Rototec-Spintec GmbH, Griesheim, Germany). Prior to inserting the vessel into the resonance cavity of the EPR spectrometer, feathers were pressed down the quartz finger to a height of approximately 0.5 cm to ensure comparable volumes of each sample. Measurements were performed at room temperature, at X-band (9.26–9.27 GHz frequency), using the following parameters: magnetic field range 3240–3340 Gs, microwave power 1 mW, modulation frequency 100 kHz, modulation amplitude and time constant—5 Gs and 0.3 s for quantitative analysis, 1 Gs and 0.1 s for qualitative analysis. An EPR signal was recorded as its first derivative, averaged from three consecutive scans, lasting 160 s each (giving a total of 480 s of scan time per EPR spectrum). Then, the following parameters were measured: peak-to-peak amplitude, area under the microwave absorption curve (the integral intensity of the recorded signal) and linewidth of the EPR absorption curve (ΔH;59).Statistical analysesStatistical analysis was performed in R (version 4.0.2,60) using a two-way ANOVA test, with bird’s diet (control vs. PT-reduced) and applied immune challenges (LPS vs. saline-injections) as the independent variables. The following parameters were used as the dependent variables: feathers reflectance (brightness), feather growth rate, feather density (number of barbs per mm), and melanisation level (expressed as the EPR spectrum amplitude measured in arbitrary units [a.u.]). The density of melanosomes was analysed by fitting a linear mixed-effects model. In this model, melanosome density was used as the dependent variable, with diet, immunological challenge, and slice ID as independent variables, and individual ID as a random-effect term. Additionally, to assess the reliability of measurements, the intraclass correlation coefficient (i.e., technical repeatability) was calculated. The models’ residuals were checked for normality and homoscedasticity. Mean food consumption per individual was analysed by the Friedman test. Body mass before and after the experiment was analysed by fitting a linear mixed-effect model. Body mass was used as the dependent variable, whereas diet, immunological challenge, and time as the independent variables, and individual ID as a random-effect term. The model included the following interaction terms: time × diet, time × injection, and diet × injection, and was reduced by removing the non-significant interactions. Results are reported with appropriate statistical tests and estimates (accompanied by standard errors) signifying relevant factor contrasts (relative to the reference group, which in all analyses was diet: control; injection: LPS, body mass: before experiment).
Ethical noteAll applicable national and institutional guidelines for the care and use of animals were followed. The research was performed under permit no. 25/2019 (with a supplementary permit no. 78/2020) from the 2nd Local Institutional Animal Care and Use Committee in Kraków. More

RESEARCH BRIEFINGS
02 November 2022

The world’s largest tropical peatland complex is in the central Congo Basin. A drying of the climate between 5,000 and 2,000 years ago triggered decomposition of peat in the Congo Basin and emission of carbon into the atmosphere. The tipping point at which drought results in carbon release might accelerate future climate change if regional droughts become more common. More

Bobbink, R. et al. Ecol. Appl. 20, 30–59 (2010).Article 
PubMed 

Google Scholar 
Olff, H. & Ritchie, M. E. Trends Ecol. Evol. 13, 261–265 (1998).Article 
PubMed 

Google Scholar 
DeMalach, N., Zaady, E. & Kadmon, R. Ecol. Lett. 20, 60–69 (2017).Article 
PubMed 

Google Scholar 
Borer, E. T. et al. Nature 508, 517–520 (2014).Article 
PubMed 

Google Scholar 
Harpole, W. S. et al. Nature 537, 93–96 (2016).Article 
PubMed 

Google Scholar 
Eskelinen, A., Harpole, W. S., Jessen, M.-T., Virtanen, R. & Hautier, Y. Nature https://doi.org/10.1038/s41586-022-05383-9 (2022).Article 

Google Scholar 
Koerner, S. E. et al. Nature Ecol. Evol. 2, 1925–1932 (2018).Article 
PubMed 

Google Scholar 
Chesson, P. Annu. Rev. Ecol. Syst. 31, 343–366 (2000).Article 

Google Scholar 
Coley, P. D., Bryant, J. P. & Chapin, F. S. Science 230, 895–899 (1985).Article 
PubMed 

Google Scholar 
Hautier, Y., Niklaus, P. A. & Hector, A. Science 324, 636–638 (2009).Article 
PubMed 

Google Scholar 
Allan, E. & Crawley, M. J. Ecol. Lett. 14, 1246–1253 (2011).Article 
PubMed 

Google Scholar  More

GPP dataAs our primary productivity product we used the GOSIF GPP dataset21 which utilizes the linear relationship between GPP and remotely-sensed SIF34. GOSIF GPP is available globally at 0.05° spatial resolution for the period 2000–2021, with the period 2001–2015 selected here (for a short summary of all datasets used in this study see Supplementary Table 3). GOSIF GPP is based on a gridded SIF product (GOSIF)34 which uses MODIS enhanced vegetation index and meteorological data for spatial scaling and is trained with millions of SIF observations from the coarser-resolution Orbiting Carbon Observatory-235. The global coverage of GOSIF and the close relationship between SIF and GPP allow for an independent assessment of how land cover changes affect GPP in different regions around the world. For instance, SIF has been shown to capture the high GPP in the US Corn Belt derived from flux towers, while ecosystem models underestimated it36. While GPP can thus be empirically estimated from satellite SIF observations relatively reliably (even though some assumptions like the linear GPP–SIF relationship and its universality across biomes are still debated20,37,38,39), the calculation of NPP needs additional assumptions of autotrophic respiration. Therefore, we focused our study on GPP, but we included an NPP product in our uncertainty analysis. In addition to that, to account for the challenges and uncertainties in global GPP estimates we included four alternative GPP products in our sensitivity analysis (see below).Land cover mappingGridded land cover was derived from ESA-CCI22, a global land cover product designed for climate science. ESA-CCI is available at 300 m spatial resolution for the 1992–2020 period (https://cds.climate.copernicus.eu/). We first classified ESA-CCI land covers to forests, grasslands, and croplands according to IPCC classification: classes 50–100, 160, 170 forests (2,022,283 grid cells); classes 110 and 130 grasslands (509,297 grid cells); classes 10–40 croplands (950,025 grid cells). We focus on these three major land cover types to facilitate our analysis. We then converted the resulting map to 0.05° resolution by determining the prevalent (i.e., mode) land cover for each grid cell using the aggregate function from the raster package40 and only included grid cells in our training data in which the prevalent land cover was constant over the period 2001–2015. Other classes (e.g., cropland/natural vegetation mosaics) and grid cells where the land cover changed over the 2001–2015 period were not used for the RF training.Random forestsRF is a popular and efficient supervised machine learning technique which can be applied for classification and regression problems41. While complex, it is still easier to interpret compared to other machine learning methods such as Artificial Neural Networks. It has recently been applied to a wide range of ecological research questions, including the prediction of food42 and bioenergy43 crop yields, potential natural vegetation31, forest aboveground biomass44, soil respiration45, and soil carbon emissions from land-use change5 and is thus well suited for our approach. The “Forests” refer to a number of individual decision trees. For each tree, a random sample of the training data is selected and split multiple times based on a random subset of variables from which the one minimizing the weighted variance is selected by the algorithm. Model performance is computed directly on out-of-bag (OOB) data which is randomly omitted from the training data (36.8% of all grid cells). When RF is applied to new data, a weighted prediction of each individual decision tree contributes to the overall prediction. Variance in the individual trees, e.g., by selecting random subsets of the observations and random variables at each node improves the overall RF predictive skill. Model training and prediction were done using the R ranger package46. After initial testing (see Supplementary Fig. S11) we decided to set the number of individual decision trees to 800 and the number of variables to possibly split at in each node to 10. As the good evaluation measures of RF algorithms can be related to spatial autocorrelation24 we also tested a coordinate-only model and performed a leave-one-out cross validation including spatial buffers (Supplementary Discussion 2, Supplementary Fig. S3). Due to the large computational effort we reduced the number of decision trees to 100 for the buffered leave-one-out cross validation.Predictor variablesWe predicted forest, grassland, and cropland potential GPP using the following 20 predictor variables in our RF algorithm: mean annual surface temperature (Tmean), mean diurnal temperature range (Tdiurnal), temperature seasonality (Tseason; standard deviation), minimum temperature of the coldest month (Tmin), annual temperature range (Tannual), mean temperature of the warmest quarter (Twarmest), mean annual precipitation (Pmean), precipitation seasonality (Pseason; coefficient of variation), precipitation of the wettest quarter (Pwettest), precipitation of the driest quarter (Pdriest), precipitation of the warmest quarter (Pwarmest), mean annual solar radiation (SR), growing degree days (GDD), relative humidity (RH), soil clay content (Clay), elevation (EL), nitrogen deposition (Ndep), nitrogen fertilization (NF), pesticide application (Pest), and gross domestic product (GDP; a proxy for agricultural management input other than NF and Pest). Overall Tmean, Tannual, and Pmean were the most important predictor variables (see Supplementary Discussion 3 and Fig. S12). We also tested other predictors (including additional bioclimatic variables, soil pH, irrigation, or phosphate fertilization) but found only negligible improvements in RF evaluation metrics and hence decided to restrict our analysis to the 20 predictors mentioned above.Climate variables were taken from the CHELSA dataset47,48, remapped to 0.05° spatial resolution using the aggregate function from the raster package40. To only include years overlapping with our GPP data we used the CHELSA time-series data for the 2001–2013 period if available and 1979–2013 climatologies elsewise. Clay was derived from the Regridded Harmonized World Soil Database v1.249. Ndep was taken from ISIMIP2b50, bilinear remapped from 0.5° to 0.05° spatial resolution using Climate Data Operators33. Elevation was obtained from WorldClim51. NF and Pest were derived from country-specific FAO data (e.g., https://ourworldindata.org/grapher/pesticide-use-per-hectare-of-cropland), i.e., we used the same value for all grid cells in a country. GDP was obtained from ref.52.Suitable areaFor the comparison of potential forest, grassland, and cropland GPP in Fig. 1g–i we only included grid cells suitable for all three land cover types. For forests, we assumed forest cover possible if the grid cell is currently forested (e.g., all grid cells of our forest training data) or if the potential natural forest cover exceeds 36.3%. This threshold represents the 5th percentile of all currently forested grid cells. Potential natural forest cover was derived from a potential natural vegetation map, available for 20 biomes at 0.00833° spatial resolution31. To convert these biomes into potential natural forest cover we assumed 100% forest cover for the ten forest biomes and 30% forest cover for tropical savannah. Other biomes were not considered. We then aggregated the map to 0.05° spatial resolution by computing the mean of 36 grid cells using the aggregate function form the raster package40 (see Supplementary Fig. S5 for the resulting map). For grasslands and croplands, we computed the 5th percentile of Tmean and Pmean in the training data (− 9.9 °C and 165 mm for grasslands and 2.7 °C and 295 mm for croplands, respectively) and removed all grid cells below those thresholds, assuming these areas to be too cold or too dry for the respective land cover type. Finally, we calculated the land cover with the highest potential GPP for all overlapping grid cells.Sensitivity analysisTo explore the sensitivity and uncertainty of our RF approach we repeated our prediction using different input datasets, potential forest cover, and machine-learning approaches. The importance of the underlying potential forest map was estimated by replacing our potential forest map (Supplementary Fig. S5) by the LUH2 potential forest map (Supplementary Fig. S13)23. To explore the dependency on the land cover product we repeated our RF prediction using the spatially aggregated MODIS land cover map (MCD12C1; IGBP scheme), available at 0.05° spatial resolution53. We classified grid cells of classes 1, 2, 3, 4, 5, (all forests), 8 (woody savannahs) and 9 (savannahs) as forest. Classes 8 and 9 were included in forest because otherwise forest cover would be underestimated in the temperate and boreal zone. Class 10 was classified as grassland and class 12 as cropland. A comparison of ESA-CCI with MODIS reveals a substantially larger cropland area in ECA-CCI but a smaller grassland area (Supplementary Fig. S14).The sensitivity to the climate product was tested by repeating our analysis using predictor variables from the WorldClim climatologies (1970–2000)51, aggregated from 30 s to 0.05° spatial resolution using the aggregate function from the raster package40. In contrast to CHELSA, growing degree days and relative humidity were not available from WorldClim but we included water vapour pressure as additional predictor.We also tested four alternative global GPP products. The vegetation photosynthesis model (VPM) product, available for the period of interest at 0.05° spatial resolution, is based on improved light use efficiency theory and is driven by remotely sensed datasets and reanalysis climate data and land cover classification which also distinguishes C3 vs. C4 photosynthesis pathways54. The second product is derived from remote sensing considering radiation and canopy conductance limitations on GPP and is available at 0.05° resolution for the 2001–2012 period55. Land cover is not an input variable. The third product, FLUXCOM, uses machine learning to scale FLUXNET site GPP to the globe56,57. FLUXCOM is available at 0.0833° resolution and was conservative remapped to 0.05° using Climate Data Operators33 meaning that the GPP of different land cover types might be mixed in regions with heterogeneous land cover patterns. The forth product is the MODIS MOD17A3 GPP product58, available for the 2001–2013 period and aggregated to 0.05° resolution using the raster package40. It is derived from meteorological data, fraction of absorbed photosynthetic active radiation/leaf area index, and land cover. As there is also a MOD17A3 NPP product available we additionally conducted a prediction for potential NPP. The MOD17A3 NPP product is calculated as GPP minus maintenance and growth respiration estimated from allometric relationships linking daily biomass and annual growth of plant tissues to leaf area index58. This leads to additional uncertainty compared to the MOD17A3 GPP product.To test the effect of an alternative RF algorithm we repeated our prediction with the RF algorithm from the Python scikit-learn library59 using the same number of decision trees (800). Additionally, we tested another machine-learning technique, a deep neural network (DNN), using the PyTorch library60. We selected 10 linear layers with 5 times alternating 128 and 256 nodes and a sigmoid output function. All layers were connected using the rectified linear unit activation function. We used the adamW optimizer with 0.0003 learning rate and 2000 epochs of training. To prevent overfitting, we included a 10% dropout after the 7th layer. Lastly, we included a very simple estimate of the most productive land cover based on the nearest neighbour using scikit-learn’s BallTree implementation together with the Haversine formula. For each grid cell we searched for the nearest forest, grassland, and cropland grid cell and assigned the respective GPP also to this grid cell. We thus assumed that environmental conditions are more or less identical in these grid cells, which might be a reasonable assumption for many locations but less reliable in complex terrain or in large homogeneous regions like the central Amazon rainforest where the nearest cropland/grassland grid cell might be located far away.Land-use change scenariosTo estimate the effects of historical and potential future land cover changes on global GPP we applied LUH2 scenarios23 which also serve as input data for ESMs participating in CMIP6. Land-use changes over the historical period are based on the HYDE reconstruction3, while future projections were developed by different Integrated Assessment Models combining various assumptions of socio-economic behaviour (SSPs) with climate mitigation targets (RCPs). Annual fractions for the two land cover classes cropland (sum of 5 crop types) and managed grassland (sum of pasture and rangeland) were available for each scenario at 0.25° resolution (https://luh.umd.edu/). We converted to 0.05° resolution assuming the same land cover fractions for all 25 grid cells around the LUH2 grid cells. We considered the following land cover transitions: forest to managed grassland, forest to cropland, and natural grassland to cropland (and reverse transitions for future scenarios). Transitions in areas suitable for only two land cover types were also included. Conversions of natural grasslands to managed grasslands were assumed not to affect productivity. We assumed the original land cover of a grid cell to be either forest (i.e., potential forest cover  > 36.3%) or natural grassland and accordingly multiplied the converted areas by the differences in potential GPP derived from our RF approach. Our broad forest definition including open tree cover (see above) and the fact that we assumed a change from 100 to 0% forest area in deforested grid cells results in a total historical deforestation area substantially larger than estimated in a recent study (2.4 Mkm2 vs. 1.6 Mkm2)61. These assumptions, however, do not impair our GPP estimate as our approach implicitly accounts for gradients in forest productivity (open forests tend to have lower GPP than closed forests). To test the sensitivity of the resulting GPP reduction we also applied the potential GPP maps from our uncertainty analysis to historical land-use changes (Supplementary Fig. S6). For future land cover changes we investigated changes over the 2015–2100 period for all available LUH2 scenarios: SSP1-1.9, SSP2-2.6, SSP4-3.4, SSP5-3.4, SSP2-4.5, SSP4-6.0, SSP3-7.0, and SSP5-8.5. Land-use activities in these scenarios range from large-scale deforestation (e.g., SSP3-7.0) to reforestation (e.g., SSP1-1.9) (Supplementary Fig. S7).Earth System ModelsWe compared the potential GPP estimated by our RF algorithm to simulations of eight ESMs participating in CMIP6 (CESM2-CLM562, CNRM-ESM2.1-Surfex 8.0c63, EC-Earth3-Veg-LPJ-GUESSv464, GFDL-ESM4-GFDL-LM4.165, IPSL-CM6A-LR-ORCHIDEEv2.066, MIROC-ES2L-MATSIRO6.0 + VISIT-e ver.1.067, MPI-ESM1-2-LR-JSBACH3.2068, UKESM1-0-LL-JULES-ES-1.069) with an explicit representation of natural vegetation and at least one agricultural land cover class (cropland or managed grassland) in their vegetation sub-model. We selected these ESMs so that all vegetation models implemented in more than one ESM were represented only once (e.g., the JSBACH vegetation model is a component of both MPI-ESM1-2-LR and AWI-ESM). For each ESM, the variable gppLut was downloaded from the CMIP6 archive (https://esgf-data.dkrz.de/search/cmip6-dkrz/) for the historical simulations. These files contain simulated GPP for natural vegetation, pasture, and cropland for which we calculated the 2001–2014 mean (2014 is the last year of the historical period). ESMs use fractional land covers for each grid cell, meaning that climatic drivers are inherently the same for all land cover types within a grid cell and simulated productivities can therefore be directly compared. As ESMs differ in their spatial resolution we bilinear remapped all output to 0.05° resolution using Climate Data Operators33 to allow for a fair comparison across models. To assess the sensitivity of our results to the interpolation method we also tested conservative remapping which results in slightly different maps (Supplementary Fig. S15) and usually larger model biases (Supplementary Table 2). In addition, ESMs differ in where they simulate forests in natural vegetation areas, and therefore we removed all grid cells from the comparison where at least one ESM simulated no tree productivity/cover/biomass in order to avoid comparing the GPP of natural grasslands to managed grasslands. We provide maps based on the original output for each ESM in Supplementary Fig. S10.FLUXNET dataWe compared our predictions of potential GPP to FLUXNET Tier 1 eddy covariance measurements (Supplementary Fig. S16)70. We included all forest, woody savannah (classified as forest), grassland and cropland sites21 which were located in suitable areas for the respective land cover. Mean GPP was calculated as the mean of the GPP estimates based on the night-time (GPP_NT_VUT_REF) and day-time (GPP_DT_VUT_REF) partitioning method. As some sites only had a few years of data, all available years were considered (i.e., site mean GPP was calculated for a different time period than 2001–2015). Comparisons were made with the potential GPP in the respective grid cell in which the site was located (i.e., not calibrated to site conditions). More

COVID variant family expandsSince the Omicron variant of SARS-CoV-2 emerged in late 2021, it has spawned a series of subvariants that have sparked global waves of infection. In the past few months, scientists have identified more than a dozen extra subvariants to watch. There are so many that they’re being called a swarm, or ‘variant soup’. BQ.1.1 (a descendant of BQ.1) and XBB seem to be rising to the top, possibly because they have many mutations in a key region of the viral spike protein called the receptor binding domain, which is required to infect cells.

Source: NextStrain

The variants near youIn Europe and North America, SARS-CoV-2 variants in the BQ.1 family are rising quickly and are likely to drive infection waves as these regions enter winter. They are also a common ingredient of the variant soup in South Africa, Nigeria and elsewhere in Africa. XBB, by contrast, looks likely to dominate infections in Asia; it recently drove a wave of infections in Singapore.

Source: Moritz Gerstung, Cov-Spectrum.org and GISAID

Money worries for science studentsEighty-five per cent of graduate students who responded to a Nature survey are worried about the increasing cost of living, and 25% are concerned about their growing student debt. Forty-five per cent said that rising inflation could cause them to reconsider whether to continue their science studies. The survey involved more than 3,200 self-selected PhD and master’s students from around the world.

How species suffer in heatwavesEven a small temperature rise has a severe effect on animal mortality, and understanding this relationship is important for predicting the effects of heatwaves caused by climate change. A paper in Nature used published data to examine how changes in temperature affect the rate of biological processes, such as movement or metabolism, at permissive temperatures — those at which species function normally. They also looked at how higher, stressful temperatures affect the rate of heat failure (irreversible heat injuries that result in death). This graph shows that rising temperatures drive a very rapid increase in heat-failure rates in frogs and molluscs. These high sensitivities suggest that when there is no way to escape hot conditions, species can quickly succumb. More

Woolhouse MEJ, Gowtage-Sequeria S. Host range and emerging and reemerging pathogens. Emerg Infect Dis. 2005;11:1842–7.Article 
PubMed 
PubMed Central 

Google Scholar 
Allen T, Murray KA, Zambrana-Torrelio C, Morse SS, Rondinini C, Di Marco M, et al. Global hotspots and correlates of emerging zoonotic diseases. Nat Commun. 2017;8:1–10.Article 
CAS 

Google Scholar 
Van Kerkhove MD, Ly S, Holl D, Guitian J, Mangtani P, Ghani AC, et al. Frequency and patterns of contact with domestic poultry and potential risk of H5N1 transmission to humans living in rural Cambodia. Influenza Other Respir Viruses. 2008;2:155–63.Article 
PubMed 
PubMed Central 

Google Scholar 
Gaythorpe KAM, Hamlet A, Cibrelus L, Garske T, Ferguson NM. The effect of climate change on yellow fever disease burden in Africa. eLife. 2020;9:1–27.Article 

Google Scholar 
Faust CL, McCallum HI, Bloomfield LSP, Gottdenker NL, Gillespie TR, Torney CJ, et al. Pathogen spillover during land conversion. Ecol Lett. 2018;21:471–83.Article 
PubMed 

Google Scholar 
Gog J, Woodroffe R, Swinton J. Disease in endangered metapopulations: The importance of alternative hosts. Proc R Soc B Biol Sci. 2002;269:671–6.Article 

Google Scholar 
Jones BA, Grace D, Kock R, Alonso S, Rushton J, Said MY, et al. Zoonosis emergence linked to agricultural intensification and environmental change. Proc Natl Acad Sci USA. 2013;110:8399–404.Article 
CAS 
PubMed 
PubMed Central 

Google Scholar 
White LA, Forester JD, Craft ME. Disease outbreak thresholds emerge from interactions between movement behavior, landscape structure, and epidemiology. Proc Natl Acad Sci USA. 2018;115:7374–9.Article 
CAS 
PubMed 
PubMed Central 

Google Scholar 
Altizer S, Bartel R, Han BA. Animal migration and infectious disease risk. Science. 2011;331:296–302.Article 
CAS 
PubMed 

Google Scholar 
Ludwig SC, Roos S, Bubb D, Baines D. Long-term trends in abundance and breeding success of red grouse and hen harriers in relation to changing management of a Scottish grouse moor. Wildl Biol. 2017;2017:wlb.00246.Article 

Google Scholar 
Newton I. Weather-related mass-mortality events in migrants. Ibis. 2007;149:453–67.Article 

Google Scholar 
Ropert-Coudert Y, Kato A, Meyer X, Pellé M, MacIntosh AJJ, Angelier F, et al. A complete breeding failure in an Adélie penguin colony correlates with unusual and extreme environmental events. Ecography. 2015;38:111–3.Article 

Google Scholar 
Newmark WD, Stanley TR. Habitat fragmentation reduces nest survival in an Afrotropical bird community in a biodiversity hotspot. Proc Natl Acad Sci USA. 2011;108:11488–93.Article 
CAS 
PubMed 
PubMed Central 

Google Scholar 
Tuyttens FaM, Macdonald DW, Rogers LM, Cheeseman CL, Roddam AW. Comparative study on the consequences of culling badgers (Meles meles) on biometrics, population dynamics and movement. J Anim Ecol. 2000;69:567–80.Article 

Google Scholar 
Frafjord K. Influence of reproductive status: Home range size in water voles (Arvicola amphibius). PLoS ONE. 2016;11:1–13.Article 

Google Scholar 
Begg CM, Begg KS, Du Toit JT, Mills MGL. Spatial organization of the honey badger Mellivora capensis in the southern Kalahari: Home-range size and movement patterns. J Zool. 2005;265:23–35.Article 

Google Scholar 
Bronikowski AM, Cords M, Alberts SC, Altmann J, Brockman DK, Fedigan LM, et al. Female and male life tables for seven wild primate species. Sci Data. 2016;3:1–8.Article 

Google Scholar 
Mitchell GW, Woodworth BK, Taylor PD, Norris DR. Automated telemetry reveals age specific differences in flight duration and speed are driven by wind conditions in a migratory songbird. Mov Ecol. 2015;3:1–13.Article 

Google Scholar 
Frankish CK, Manica A, Phillips RA. Effects of age on foraging behavior in two closely related albatross species. Mov Ecol. 2020;8:1–17.Article 

Google Scholar 
Tirpak JM, Giuliano WM, Allen TJ, Bittner S, Edwards JW, Friedhof S, et al. Ruffed grouse-habitat preference in the central and southern Appalachians. Ecol Manag. 2010;260:1525–38.Article 

Google Scholar 
Zhu WW, Garber PA, Bezanson M, Qi XG, Li BG. Age- and sex-based patterns of positional behavior and substrate utilization in the golden snub-nosed monkey (Rhinopithecus roxellana). Am J Primatol. 2015;77:98–108.Article 
PubMed 

Google Scholar 
Tian H, Yu P, Bjørnstad ON, Cazelles B, Yang J, Tan H, et al. Anthropogenically driven environmental changes shift the ecological dynamics of hemorrhagic fever with renal syndrome. PLOS Pathog. 2017;13:e1006198.Article 
PubMed 
PubMed Central 

Google Scholar 
George DB, Webb CT, Farnsworth ML, O’Shea TJ, Bowen RA, Smith DL, et al. Host and viral ecology determine bat rabies seasonality and maintenance. Proc Natl Acad Sci USA. 2011;108:10208–13.Article 
CAS 
PubMed 
PubMed Central 

Google Scholar 
van Dijk JG, Verhagen JH, Wille M, Waldenström J. Host and virus ecology as determinants of influenza A virus transmission in wild birds. Curr Opin Virol. 2018;28:26–36.Article 
PubMed 

Google Scholar 
Chong R, Shi M, Grueber CE, Holmes EC, Hogg CJ, Belov K, et al. Fecal Viral Diversity of Captive and Wild Tasmanian Devils Characterized Using Virion-Enriched Metagenomics and Metatranscriptomics. J Virol. 2019;93:e00205–19.Article 
CAS 
PubMed 
PubMed Central 

Google Scholar 
François S, Pybus OG. Towards an understanding of the avian virome. J Gen Virol. 2020;101:785–90.Article 
PubMed 
PubMed Central 

Google Scholar 
Springer A, Fichtel C, Al-Ghalith GA, Koch F, Amato KR, Clayton JB, et al. Patterns of seasonality and group membership characterize the gut microbiota in a longitudinal study of wild Verreaux’s sifakas (Propithecus verreauxi). Ecol Evol. 2017;7:5732–45.Article 
PubMed 
PubMed Central 

Google Scholar 
Aivelo T, Laakkonen J, Jernvall J. Population-and individual-level dynamics of the intestinal microbiota of a small primate. Appl Environ Microbiol. 2016;82:3537–45.Article 
CAS 
PubMed 
PubMed Central 

Google Scholar 
van Dongen WFD, White J, Brandl HB, Moodley Y, Merkling T, Leclaire S, et al. Age-related differences in the cloacal microbiota of a wild bird species. BMC Ecol. 2013;13:11.Cleaveland S, Laurenson MK, Taylor LH. Diseases of humans and their domestic mammals: Pathogen characteristics, host range and the risk of emergence. Philos Trans R Soc B Biol Sci. 2001;356:991–9.Article 
CAS 

Google Scholar 
Wille M, Shi M, Hurt AC, Klaassen M, Holmes EC. RNA virome abundance and diversity is associated with host age in a bird species. Virology. 2021;561:98–106.Article 
CAS 
PubMed 

Google Scholar 
Negrey JD, Thompson ME, Langergraber KE, Machanda ZP, Mitani JC, Muller MN, et al. Demography, life-history trade-offs, and the gastrointestinal virome of wild chimpanzees. Philos Trans R Soc B Biol Sci. 2020;375:20190613.Article 

Google Scholar 
Hill SC, Manvell RJ, Schulenburg B, Shell W, Wikramaratna PS, Perrins C, et al. Antibody responses to avian influenza viruses in wild birds broaden with age. Proc R Soc B Biol Sci. 2016;283:20162159.Article 

Google Scholar 
Perrins CM, Ogilvie MA. A study of the Abbotsbury mute swans (Cygnus olor). Wildfowl. 1981;32:35–47.
Google Scholar 
Perrins CM, McCleery RH, Ogilvie MA. A study of the breeding Mute Swans Cygnus olor at Abbotsbury. Wildfowl. 1994;45:1–14.
Google Scholar 
Perrins C. Survival rates of young mute swans Cygnus olor. Wildfowl Suppl. 1991;45:95–103.
Google Scholar 
McCleery RH, Perrins C, Wheeler D, Groves S. Population structure, survival rates and productivity of mute swans breeding in a colony at Abbotsbury, Dorset, England. Waterbirds Waterbird Soc. 2002;25:201.
Google Scholar 
Matrozis R A 30-year (1988–2017) study of Mute Swans Cygnus olor in Riga, Latvia. Wildfowl. 2019;14:164–77.Charmantier A, Perrins C, McCleery RH, Sheldon BC. Quantitative genetics of age at reproduction in wild swans: Support for antagonistic pleiotropy models of senescence. Proc Natl Acad Sci USA. 2006;103:6587–92.Article 
CAS 
PubMed 
PubMed Central 

Google Scholar 
Hill SC, Hansen R, Watson S, Coward V, Russell C, Cooper J, et al. Comparative micro-epidemiology of pathogenic avian influenza virus outbreaks in a wild bird population. Philos Trans R Soc B Biol Sci. 2019;374:20180259.Cotten M, Oude Munnink B, Canuti M, Deijs M, Watson SJ, Kellam P, et al. Full genome virus detection in fecal samples using sensitive nucleic acid preparation, deep sequencing, and a novel iterative sequence classification algorithm. PLoS ONE. 2014;9:e93269.Article 
PubMed 
PubMed Central 

Google Scholar 
Boom R, Sol CJA, Salimans MMM, Janses CL, Wertheim Van Dillen PME, Van Der Noordaa J. Rapid and simple method for purification of nucleic acids R. J Clin Microbiol. 1990;28:495–503.Article 
CAS 
PubMed 
PubMed Central 

Google Scholar 
Endoh D, Mizutani T, Kirisawa R, Maki Y, Saito H, Kon Y, et al. Species-independent detection of RNA virus by representational difference analysis using non-ribosomal hexanucleotides for reverse transcription. Nucleic Acids Res. 2005;33:1–11.Article 

Google Scholar 
Martin M. Cutadapt removes adapter sequences from high-throughput sequencing reads. EMB. 2011;17:10–12.
Google Scholar 
Zerbino DR, Birney E. Velvet: algorithms for de novo short read assembly using de Bruijn graphs. Genome Res. 2008;18:821–9.Article 
CAS 
PubMed 
PubMed Central 

Google Scholar 
Buchfink B, Xie C, Huson DH. Fast and sensitive protein alignment using DIAMOND. Nat Methods. 2014;12:59–60.Article 
PubMed 

Google Scholar 
Altschul SF, Gish W, Miller W, Myers EW, Lipman DJ. Basic local alignment search tool. J Mol Biol. 1990;215:403–10.Article 
CAS 
PubMed 

Google Scholar 
Langmead B Aligning short sequencing reads with Bowtie. Curr Protoc Bioinforma. 2010; Chapter 11: Unit 11.7.Kearse M, Moir R, Wilson A, Stones-Havas S, Cheung M, Sturrock S, et al. Geneious Basic: an integrated and extendable desktop software platform for the organization and analysis of sequence data. Bioinforma Oxf Engl. 2012;28:1647–9.Article 

Google Scholar 
Katoh K, Misawa K, Kuma K, Miyata T. MAFFT: a novel method for rapid multiple sequence alignment based on fast Fourier transform. Nucleic Acids Res. 2002;30:3059–66.Article 
CAS 
PubMed 
PubMed Central 

Google Scholar 
Katoh K, Standley DM. MAFFT multiple sequence alignment software version 7: Improvements in performance and usability. Mol Biol Evol. 2013;30:772–80.Article 
CAS 
PubMed 
PubMed Central 

Google Scholar 
Edgar RC. MUSCLE: multiple sequence alignment with high accuracy and high throughput. Nucleic Acids Res. 2004;32:1792–7.Article 
CAS 
PubMed 
PubMed Central 

Google Scholar 
Muhire BM, Varsani A, Martin DP SDT: A virus classification tool based on pairwise sequence alignment and identity calculation. PLoS ONE. 2014;9:e108277.Darriba D, Taboada GL, Doallo R, Posada D. ProtTest 3: fast selection of best-fit models of protein evolution. Bioinforma Oxf Engl. 2011;27:1164–5.Article 
CAS 

Google Scholar 
Stamatakis A. RAxML version 8: A tool for phylogenetic analysis and post-analysis of large phylogenies. Bioinformatics. 2014;30:1312–3.Article 
CAS 
PubMed 
PubMed Central 

Google Scholar 
Kapoor A, Simmonds P, Lipkin WI, Zaidi S, Delwart E. Use of nucleotide composition analysis to infer hosts for three novel picorna-like viruses. J Virol. 2010;84:10322–8.Article 
CAS 
PubMed 
PubMed Central 

Google Scholar 
Wood DE, Salzberg SL Kraken: Ultrafast metagenomic sequence classification using exact alignments. Genome Biol. 2014;15:R46.Lu J, Breitwieser FP, Thielen P, Salzberg SL. Bracken: estimating species abundance in metagenomics data. PeerJ Comput Sci. 2017;3:e104.Article 

Google Scholar 
Cole JR, Wang Q, Fish JA, Chai B, McGarrell DM, Sun Y, et al. Ribosomal Database Project: Data and tools for high throughput rRNA analysis. Nucleic Acids Res. 2014;42:633–42.Article 

Google Scholar 
R Core Team. R: A Language and Environment for Statistical Computing. 2019. R Foundation for Statistical Computing, Vienna, Austria.RStudio Team. RStudio: Integrated Development for R. 2015. Boston, MA, USA.Dixon P. VEGAN, a package of R functions for community ecology. J Veg Sci. 2003;14:927–30.Article 

Google Scholar 
Anderson MJ. A new method for non-parametric multivariate analysis of variance. Austral Ecol. 2001;26:32–46.
Google Scholar 
Benjamini Y, Hochberg Y. Controlling the false discovery rate: a practical and powerful approach to multiple testing. J R Stat Soc Ser B Methodol. 1995;57:289–300.
Google Scholar 
Benjamini Y, Yekutieli D. The control of the false discovery rate in multiple testing under dependency. Ann Stat. 2001;29:1165–88.Article 

Google Scholar 
Oakley BB, Lillehoj HS, Kogut MH, Kim WK, Maurer JJ, Pedroso A, et al. The chicken gastrointestinal microbiome. FEMS Microbiol Lett. 2014;360:100–12.Article 
CAS 
PubMed 

Google Scholar 
Waite DW, Taylor MW. Characterizing the avian gut microbiota: Membership, driving influences, and potential function. Front Microbiol. 2014;5:1–12.Article 

Google Scholar 
Waite DW, Taylor MW. Exploring the avian gut microbiota: Current trends and future directions. Front Microbiol. 2015;6:1–12.Article 

Google Scholar 
Zell R, Delwart E, Gorbalenya AE, Hovi T, King AMQ, Knowles NJ, et al. ICTV virus taxonomy profile: Picornaviridae. J Gen Virol. 2017;98:2421–2.Article 
CAS 
PubMed 
PubMed Central 

Google Scholar 
Cotmore SF, Agbandje-McKenna M, Canuti M, Chiorini JA, Eis-Hubinger AM, Hughes J, et al. ICTV virus taxonomy profile: Parvoviridae. J Gen Virol. 2019;100:367–8.Article 
CAS 
PubMed 
PubMed Central 

Google Scholar 
Bosch A, Guix S, Krishna N, Méndez E, Monroe SS, Pantin-Jackwood M, et al. Astroviridae. In: King A, Adams M, Carstens E, Lefkowitz E (eds). Virus taxonomy. Classification and nomenclature of viruses: ninth report of the International Committee on the Taxonomy of Viruses. 2011. Elsevier, London, pp 953-9.Risely A. Applying the core microbiome to understand host–microbe systems. J Anim Ecol. 2020;89:1549–58.Article 
PubMed 

Google Scholar 
Piepenbring AK, Enderlein D, Herzog S, Kaleta EF, Heffels-Redmann U, Ressmeyer S, et al. Pathogenesis of avian bornavirus in experimentally infected Cockatiels. Emerg Infect Dis. 2012;18:234–41.Article 
PubMed 
PubMed Central 

Google Scholar 
Anzil AP, Blinzinger K, Mayr A. Persistent Borna virus infection in adult hamsters. Arch Für Gesamt Virusforsch. 1973;40:52–57.Article 
CAS 

Google Scholar 
Heffels-Redmann U, Enderlein D, Herzog S, Piepenbring A, Bürkle M, Neumann D, et al. Follow-Up Investigations on Different Courses of Natural Avian Bornavirus Infections in Psittacines. Avian Dis. 2012;56:153–9.Article 
PubMed 

Google Scholar 
Rubbenstroth D, Brosinski K, Rinder M, Olbert M, Kaspers B, Korbel R, et al. No contact transmission of avian bornavirus in experimentally infected cockatiels (Nymphicus hollandicus) and domestic canaries (Serinus canaria forma domestica). Vet Microbiol. 2014;172:146–56.Article 
PubMed 

Google Scholar 
Olsen I The Family Fusobacteriaceae. In: Rosenberg E, Delong EF, Lory S, Stackebrandt E, Thompson F (eds). The Prokaryotes: Firmicutes and Tenericutes, 4th ed. 2014. pp 109-32.Imhoff JF The Family Chlorobiaceae. In: Rosenberg E, Delong EF, Lory S, Stackebrandt E, Thompson F (eds). The Prokaryotes: Other Major Lineages of Bacteria and The Archaea, 4th ed. 2014. pp 501-14.Cho JC The Family Lentisphaeraceae. In: Rosenberg E, Delong EF, Lory S, Stackebrandt E, Thompson F (eds). The Prokaryotes: Other Major Lineages of Bacteria and The Archaea, 4th ed. 2014. pp 705-10.Karami A, Sarshar M, Ranjbar R, Zanjani RS The Phylum Spirochaetaceae. In: Rosenberg E, Delong EF, Lory S, Stackebrandt E, Thompson F (eds). The Prokaryotes: Other Major Lineages of Bacteria and The Archaea, 4th ed. 2014. pp 915-29.McBride MJ The Family Flavobacteriaceae. In: Rosenberg E, Delong EF, Lory S, Stackebrandt E, Thompson F (eds). The Prokaryotes: Other Major Lineages of Bacteria and The Archaea, 4th ed. 2014. pp 643-76.Van Dijk JGB, Hoye BJ, Verhagen JH, Nolet BA, Fouchier RAM, Klaassen M. Juveniles and migrants as drivers for seasonal epizootics of avian influenza virus. J Anim Ecol. 2014;83:266–75.Article 
PubMed 

Google Scholar 
Chevalier V, Marsot M, Molia S, Rasamoelina H, Rakotondravao R, Pedrono M, et al. Serological evidence of West Nile and Usutu viruses circulation in domestic and wild birds in wetlands of Mali and Madagascar in 2008. Int J Environ Res Public Health. 2020;17:1998.Guy, JS Turkey Viral Hepatitis. Diseases of Poultry, 12th Edition. 2008. Wiley Blackwell, pp 426-30.Davies ZG, Fuller RA, Dallimer M, Loram A, Gaston KJ. Household factors influencing participation in bird feeding activity: a national scale analysis. PLOS ONE. 2012;7:e39692.Article 
CAS 
PubMed 
PubMed Central 

Google Scholar 
Shutt JD, Trivedi UH, Nicholls JA. Faecal metabarcoding reveals pervasive long-distance impacts of garden bird feeding. Proc R Soc B Biol Sci. 2021;288:20210480.Article 

Google Scholar 
Minich JJ, Sanders JG, Amir A, Humphrey G, Gilbert JA, Knight R. Quantifying and understanding well-to-well contamination in microbiome research. mSystems. 2019;4:e00186–19.Article 
PubMed 
PubMed Central 

Google Scholar  More