Philippa Kaur

Soil organic matter is essential for maintaining soil health and sustaining plant growth. Loss of soil organic matter often leads to degradation of soil quality1. It also constitutes the largest terrestrial organic carbon (C) pool (~ 2,400 Pg C in top 2 m soil). This soil organic C (SOC) pool is three times greater than the amount of C in the atmosphere2. An increase or decrease of C in soils by only a small percent represents a substantial C sink or source for atmospheric CO2. Studies have been conducted to predict SOC of croplands under climate change3,4,5. However, crop management changes with adaptation to future climate6,7,8,9 (such as changes in varieties and planting dates) were ignored in most studies. Thus, we carried out a regional study to assess the impact of projected climate change and elevated CO2 on SOC in agricultural systems with management adaptation.
A change in SOC is a result of the net effect of the changes in SOC decomposition rates (the main C outflow) and C inputs from plants10 (main C inflow). In a warmer climate, the higher temperature could increase the decomposition rate in both the short and long term11. Decomposition will also be sensitive to increases or decreases in precipitation10 predicted by climate models8. Carbon input is mainly from root and surface residue litter that is not removed from the system through harvest, grazing or burning. Changes in temperature, precipitation, and atmospheric CO2 all affect this input through plant growth and production6. These climatic changes are also likely to affect management practices in the future7. Although systemic management changes are possible that can limit the impact of climate change, such as moving from dryland to irrigated systems, smaller adjustments including changes in varieties and planting dates have less barriers for adoption. These adjustments to management have been found to affect both crop production (C input) and decomposition7,12, with potentially larger effects on crop production7. These interactions are complex, and the overall change to SOC pools in agricultural lands remains uncertain.
Here we present the results from a simulation of SOC dynamics for more than 54,000 locations, covering 34,600,000 ha of cropland in the U.S. Corn Belt (Supplementary Fig. 1a). These locations have historically been managed with a corn (Zea mays L.) and soybean (Glycine max L.) rotation, which is the most common crop rotation in this region. Corn and soybean yields in this region account for about 85% of US crop production13. In our simulation, the widely used biogeochemical model DayCent14,15 was driven by weather data from the Representative Concentration Pathway 4.5 (RCP 4.5) climate change scenario8. Predictions from three General Circulation Models (GCMs; GFDL-CM3, MIROC-ESM, and MRI-CGCM3) were downscaled to generate daily weather (32 km grid) and used to assess the uncertainty. The main goal was to quantify the change in surface soil C (0–20 cm) under management adaptation with the selection of alternative crop varieties and planting dates to maintain high levels of crop production in future climate scenarios from 2041 to 2071. Alternative crop varieties were based on those currently available in the United States without consideration of additional crop breeding or genetic modifications that could further enhance production in the future. A no climate change scenario (historical weather with current CO2 level) was used as the baseline for comparisons.
Crop production trends
We found that, without adaptation, the average corn yield during the 2041–2071 period in the Corn Belt would drop in all three GCM climate scenarios in comparison with the baseline with no climate change (Fig. 1). Yield decreased by 17% (GFDL-CM3), 34% (MIROC-ESM), and 2% (MRI-CGCM3) under the respective scenarios. The difference between GCMs can be attributed to differences in the projected temperature and precipitation (Supplementary Fig. 2). Similar yield decreases were seen for soybean as found for corn if there is no adaptation (decreased by 13% for GFDL-CM3, 28% for MIROC-ESM, and increased by 8% for MRI-CGCM3).
Figure 1

Predicted grain yield of (a) corn and (b) soybean in the 2041 to 2070 period.

Full size image

With adaptation, both corn and soybean yields increased compared with no adaptation: an increase of 5% for corn and 19% for soybean (average of the three GCMs) compared with the baseline. The larger increase in soybean yield was due to the C3 crop being more responsive to the CO2 fertilization than C4 crops16. The standard deviation of yields across the three GCMs in the adaptation scenarios was lower than that of the no adaptation scenarios (15% and 39% lower for corn and soybean respectively). Simulated crop yields were more stable under climate change with adaptation management compared to without adaptation. The stability in yields is because, without crop adaptation associated with the selection of alternative crop varieties, an increase in temperature shortens the growing period of the crop17,18. Each GCM predicts a longer growing season with warmer temperatures. When a longer-season maturity variety was simulated as an adaptation pathway, the extended growing period allows the crop to use the full window of optimal solar radiation. This leads to similar amounts of production from year to year between GCM climate scenarios. The spatial pattern of the yield changes was very different for the two crops (Supplementary Fig. 3) due to the difference in crop response to temperature, day length (photoperiod), and elevated CO2. Our overall predictions of crop yield change were generally in agreement with other studies9,17,19,20.
Carbon input trends
Crop yields are good indicators of total net primary production and are found to be proportional to the total C input to soils in the U.S. Corn Belt21 (i.e., crop yield to total biomass does not vary much geographically). Our simulations predicted the average input to be 3.7 Mg C ha−1 year−1 under the no climate change scenario for the 2041–2070 period (Supplementary Fig. 4). With climate change but no adaptation, all counties (an administrative subdivision of a state in the U.S.) had lower C input (GCMs ensemble mean) compared with the baseline. However, with adaptation, more than half of the counties (most counties in the northern part of the region) had higher C input than those of the baseline. Compared with the no adaptation scenario, all counties in the Corn Belt maintained higher C input with adaptation. The average C input across climate scenarios in the Corn Belt region was predicted to be 3.0 and 3.9 Mg C ha−1 (a change of − 19% and 5% from the baseline) with no adaptation and adaptation scenarios, respectively. The standard deviation of the adaptation scenarios across the three GCMs was 47% lower than without adaptation, suggesting a counteracting effect of adaptation to climate change.
Decomposition factor trend
In contrast to C input, the decomposition factor, which reflected the relative change in decomposition rate due to temperature and moisture effect (ranges from 0 to 1, with higher rates associated with values closer to 1), increased in all counties regardless of the adaptation (Supplementary Fig. 4). Larger decomposition factors can be explained by an increase in soil temperature and wetter soil conditions associated with climate change projections. Wetter soil conditions were a result of a reduced transpiration rate under high levels of atmospheric CO2 and increased precipitation projected by the three GCMs. The average decomposition factor (across three GCMs) was predicted to be 0.41, 0.49, and 0.47 for baseline, without adaptation, and with adaptation scenarios. The standard deviations were 27% lower with adaptation for the decomposition factor across the three GCMs than without adaptation. If adaptation does not occur under the GCM scenarios, the growing period of the crops was reduced due to global warming and resulted in less total transpiration and wetter soil conditions, which increased decomposition. With adaptation, the longer-season maturity variety continued to grow over a longer period and consumed more water, reducing soil moisture to lower levels, thus lowering decomposition (average annual evapotranspiration was 2.0–5.1% higher in the adaptation scenario compared with no adaptation). The soil moisture differences among GCMs were lower with crop adaptation.
Soil organic C trends
Without climate change (baseline scenario), the predicted sub-region SOC in the top soil ranged from less than 30 Mg C ha−1 to more than 80 Mg C ha−1(Fig. 2a). The highest levels of SOC were found in the western part of the region where soil clay content is high (Supplementary Fig. 1). The low levels of SOC found in Michigan and northern Indiana can be attributed to the soils with high sand content. In addition, tillage intensity varies among the counties22, contributing to differences in SOC levels among counties. With climate change, there were losses of SOC in almost all counties if there was no adaptation (Fig. 2b). This was the net result of decreased C inputs and increased decomposition rates (Supplementary Fig. 4). With adaptation, the northern part of the Corn Belt tended to gain SOC (Fig. 2c) due to increased C inputs (Supplementary Fig. 4), while the other areas tended to lose SOC. Compared with no adaptation, the adaptation scenario resulted in higher SOC stocks in all counties of the Corn Belt, which was also found in a modeling study of a site in Michigan9.
Figure 2

Predicted soil organic carbon (SOC) in the top soil (0–20 cm) averaged for the three GCMs in the U.S. Corn Belt, including (a) the baseline scenario and (b) the difference between the no adaptation scenario and baseline, and (c) the difference between the adaptation scenario and baseline. Maps were generated using the R “ggmap” package39 (Version 2.6.1; https://journal.r-project.org/archive/2013-1/kahle-wickham.pdf).

Full size image

Over the period from 2041 to 2070, the average SOC in the Corn Belt reached a new equilibrium state in the baseline scenario with no climate change (Fig. 3). Without adaptation in the climate change scenarios, SOC decreased over time. Rapid loss of C was predicted for the 2041–2050 period with the rate decreasing gradually over the 2051–2070 period. With adaptation, two GCMs (GFDL-CM3 and MRI-CGCM3) predicted similar changes over time as the baseline, while the SOC values for the other GCM (MIROC-ESM) were slightly lower than the baseline.
Figure 3

Area-weighted average soil organic carbon stocks for simulations with and without adaptation for corn/soybean rotations between 2041 and 2070. These projections are based on the same historical data and initial values for SOC pools.

Full size image

These results show that adaptation with the selection of longer-season varieties can lead to a similar SOC level as the baseline and GCM scenarios under a moderate climate change projection (RCP 4.5). Although the predicted climate was very different for the three GCMs, the SOC levels with adaptation were similar and not very different from the baseline with no climate change. In contrast, without adaptation, SOC storage levels were farther apart from each other and the baseline. This indicates a strong resistance to the effects of climate change in agricultural systems in the Corn Belt region if there is the selection for alternative varieties and planting dates that are better adapted to the changing climatic conditions.
We found the variation (CV 2.89%) of the predicted total SOC stock across the GCMs in the adaptation scenarios was much smaller than the variation (CV 5.88%) in the no adaptation scenarios. This corresponds to reduced variation in C input and decomposition factors. Because the increased decomposition rate was compensated by the higher overall C input, the SOC stock maintained similar levels as the baseline without climate change.
Limitations
In this study, we did not evaluate the possibility of new technologies that could be developed in the future and influence production, decomposition and other variables influencing SOC levels, as these changes are difficult to predict. For example, new technologies associated with crop breeding and other genetic improvements, pest control, and other developments could increase C input and result in higher SOC stocks than our predictions. Management practices such as adding cover crops, which increases C input, may also further enhance SOC stocks23. In contrast, large areas of removal of corn stover for biofuel production (not considered in our simulation) could reduce the total SOC stock24. Although most SOC is concentrated in the surface soil, subsoil C has been found to respond to warming climates and affects SOC stocks25. Future research should also address subsoil C dynamics. More

1.
Moritz, M. A., Morais, M. E., Summerell, L. A., Carlson, J. M. & Doyle, J. Wildfires, complexity, and highly optimized tolerance. Proc. Natl. Acad. Sci. USA. 102, 17912–17917 (2005).
ADS  CAS  PubMed  Google Scholar 
2.
Littell, J. S. et al. Climate and wildfire area burned in western US ecoprovinces, 1916–2003. Ecol. Appl. 19, 1003–1021 (2009).
PubMed  Google Scholar 

3.
Barbero, R. et al. Multi-scalar influence of weather and climate on very large-fires in the Eastern United States. Int. J. Climatol. 35, 2180–2186 (2015).
Google Scholar 

4.
Turco, M. et al. On the key role of droughts in the dynamics of summer fires in Mediterranean Europe. Sci. Rep. 7, 81 (2017).
ADS  PubMed  PubMed Central  Google Scholar 

5.
Ruffault, J., Moron, V., Trigo, R. M. & Curt, T. Daily synoptic conditions associated with large fire occurrence in Mediterranean France: Evidence for a wind-driven fire regime. Int. J. Climatol. 37, 524–533 (2017).
Google Scholar 

6.
Abatzoglou, J. T. & Williams, A. P. Impact of anthropogenic climate change on wildfire across western US forests. Proc. Natl. Acad. Sci. 113, 11770–11775 (2016).
ADS  CAS  PubMed  Google Scholar 

7.
Turco, M. et al. Decreasing fires in mediterranean Europe. PLoS One 11, 20 (2016).
Google Scholar 

8.
Ruffault, J., Curt, T., Martin-StPaul, N. K., Moron, V. & Trigo, R. M. Extreme wildfire events are linked to global-change-type droughts in the northern Mediterranean. Nat. Hazards Earth Syst. Sci. 18, 847–856 (2018).
ADS  Google Scholar 

9.
Turco, M. et al. Climate drivers of the 2017 devastating fires in Portugal. Sci. Rep. 9, 13886 (2019).
ADS  PubMed  PubMed Central  Google Scholar 

10.
Lagouvardos, K., Kotroni, V., Giannaros, ΤΜ & Dafis, S. Meteorological conditions conducive to the rapid spread of the deadly wildfire in eastern Attica, Greece. Bull. Am. Meteorol. Soc. https://doi.org/10.1175/bams-d-18-0231.1 (2019).
Article  Google Scholar 

11.
Turco, M. et al. Exacerbated fires in Mediterranean Europe due to anthropogenic warming projected with non-stationary climate-fire models. Nat. Commun. 9, 1–9 (2018).
Google Scholar 

12.
Dupuy, J.-L. et al. Climate change impact on future wildfire danger and activity in southern Europe: A review. Ann. For. Sci. 77, 1–49 (2020).
Google Scholar 

13.
Boer, M. M. et al. Changing weather extremes call for early warning of potential for catastrophic fire. Earth’s Future 5, 1196–1202 (2017).
ADS  Google Scholar 

14.
Fernandes, P. M., Barros, A. M. G., Pinto, A. & Santos, J. A. Characteristics and controls of extremely large wildfires in the western Mediterranean Basin. J. Geophys. Res. G Biogeosci. 121, 2141–2157 (2016).
ADS  Google Scholar 

15.
Hernandez, C., Drobinski, P. & Turquety, S. How much does weather control fire size and intensity in the Mediterranean region?. Ann. Geophys. 20, 20 (2015).
Google Scholar 

16.
Jin, Y. et al. Identification of two distinct fire regimes in Southern California: Implications for economic impact and future change. Environ. Res. Lett. 10, 94005 (2015).
Google Scholar 

17.
Ruffault, J., Moron, V., Trigo, R. M. & Curt, T. Objective identification of multiple large fire climatologies: An application to a Mediterranean ecosystem. Environ. Res. Lett. 11, 075006 (2016).
ADS  Google Scholar 

18.
Duane, A., Piqué, M., Castellnou, M. & Brotons, L. Predictive modelling of fire occurrences from different fire spread patterns in Mediterranean landscapes. Int. J. Wildl. Fire 24, 407–418 (2015).
Google Scholar 

19.
Rodrigues, M., Trigo, R. M., Vega-García, C. & Cardil, A. Identifying large fire weather typologies in the Iberian Peninsula. Agric. For. Meteorol. 280, 107789 (2020).
ADS  Google Scholar 

20.
Van Wagner, C. E. Structure of the Canadian forest fire weather index. Can. For. Serv. For. Tech. Rep. 35, 37 (1987).
Google Scholar 

21.
Ruffault, J., Martin-StPaul, N., Pimont, F. & Dupuy, J.-L.L. How well do meteorological drought indices predict live fuel moisture content (LFMC)? An assessment for wildfire research and operations in Mediterranean ecosystems. Agric. For. Meteorol. 262, 391–401 (2018).
ADS  Google Scholar 

22.
Dee, D. P. et al. The ERA-Interim reanalysis: Configuration and performance of the data assimilation system. Q. J. R. Meteorol. Soc. 137, 553–597 (2011).
ADS  Google Scholar 

23.
Abatzoglou, J. T. & Kolden, C. A. Relative importance of weather and climate on wildfire growth in interior Alaska. Int. J. Wildl. Fire 20, 479–486 (2011).
Google Scholar 

24.
Gudmundsson, L., Rego, F. C., Rocha, M. & Seneviratne, S. I. Predicting above normal wildfire activity in southern Europe as a function of meteorological drought. Environ. Res. Lett. 9, 84008 (2014).
Google Scholar 

25.
Urbieta, I. R. et al. Fire activity as a function of fire–weather seasonal severity and antecedent climate across spatial scales in southern Europe and Pacific western USA. Environ. Res. Lett. 10, 114013 (2015).
ADS  Google Scholar 

26.
Paschalidou, A. K. & Kassomenos, P. A. What are the most fire-dangerous atmospheric circulations in the Eastern-Mediterranean? Analysis of the synoptic wildfire climatology. Sci. Total Environ. 539, 536–545 (2016).
ADS  CAS  PubMed  Google Scholar 

27.
Kotlarski, S. et al. Regional climate modeling on European scales: A joint standard evaluation of the EURO-CORDEX RCM ensemble. Geosci. Model Dev. 7, 1297–1333 (2014).
ADS  Google Scholar 

28.
Pereira, M. G., Trigo, R. M., Da Camara, C. C., Pereira, J. M. C. & Leite, S. M. Synoptic patterns associated with large summer forest fires in Portugal. Agric. For. Meteorol. 129, 11–25 (2005).
ADS  Google Scholar 

29.
Cardil, A., Merenciano, D. & Molina-Terrén, D. Wildland fire typologies and extreme temperatures in NE Spain. iForest. Biogeosci. For. 009, e1–e6 (2016).
Google Scholar 

30.
Belhadj-Khedher, C., El-Melki, T., Mouillot, F. Saharan hot and dry Sirocco winds drive extreme fire events in Mediterranean Tunisia (North Africa). Atmosphere 11(6), 590 (2020).
ADS  Google Scholar 

31.
Williams, A. P. et al. Temperature as a potent driver of regional forest drought stress and tree mortality. Nat. Clim. Change 3, 292–297 (2013).
ADS  Google Scholar 

32.
Allen, C. D., Breshears, D. D. & McDowell, N. G. On underestimation of global vulnerability to tree mortality and forest die-off from hotter drought in the Anthropocene. Ecosphere 6, art129 (2015).
Google Scholar 

33.
Adams, A. & Patrick, L. Temperature response surfaces for mortality risk of tree species with future drought Item Type Article Citation Temperature response surfaces for mortality risk of tree species with future drought. Environ. Res. Lett. 12, 115014 (2017).
ADS  Google Scholar 

34.
Cochard, H. A new mechanism for tree mortality due to drought and heatwaves. bioRxiv https://doi.org/10.1101/531632 (2019).
Article  Google Scholar 

35.
Brodribb, T. J., Powers, J., Cochard, H. & Choat, B. Hanging by a thread? Forests and drought. Science 368, 261–266 (2020).
ADS  CAS  PubMed  Google Scholar 

36.
Pfleiderer, P., Schleussner, C.-F., Kornhuber, K. & Coumou, D. Summer weather becomes more persistent in a 2 °C world. Nat. Clim. Change 9, 666–671 (2019).
ADS  Google Scholar 

37.
Cramer, W. et al. Climate change and interconnected risks to sustainable development in the Mediterranean. Nat. Clim. Change 8, 972–980 (2018).
ADS  Google Scholar 

38.
Alexander, M. E. & Cruz, M. G. Assessing the effect of foliar moisture on the spread rate of crown fires. Int. J. Wildl. Fire 22, 415–427 (2013).
Google Scholar 

39.
Brotons, L., Aquilué, N., de Cáceres, M., Fortin, M. J. & Fall, A. How fire history, fire suppression practices and climate change affect wildfire regimes in mediterranean landscapes. PLoS One 8, e62392 (2013).
ADS  CAS  PubMed  PubMed Central  Google Scholar 

40.
Batllori, E., Parisien, M. A., Krawchuk, M. A. & Moritz, M. A. Climate change-induced shifts in fire for Mediterranean ecosystems. Glob. Ecol. Biogeogr. 22, 1118–1129 (2013).
Google Scholar 

41.
Balch, J. K. et al. Human-started wildfires expand the fire niche across the United States. Proc. Natl. Acad. Sci. USA 114, 2946–2951 (2017).
ADS  CAS  PubMed  Google Scholar 

42.
Parisien, M. A. et al. Fire deficit increases wildfire risk for many communities in the Canadian boreal forest. Nat. Commun. 11, 20 (2020).
Google Scholar 

43.
Moreira, F. et al. Wildfire management in Mediterranean-type regions: Paradigm change needed. Environ. Res. Lett. 15, 20 (2020).
Google Scholar 

44.
Moreira, F. et al. Landscape–wildfire interactions in southern Europe: Implications for landscape management. J. Environ. Manag. 92, 2389–2402 (2011).
Google Scholar 

45.
Ruffault, J. & Mouillot, F. How a new fire-suppression policy can abruptly reshape the fire–weather relationship. Ecosphere 6, 1–19 (2015).
Google Scholar 

46.
Pimont, F., Ruffault, J., Martin-Stpaul, N. K. & Dupuy, J. L. Why is the effect of live fuel moisture content on fire rate of spread underestimated in field experiments in shrublands?. Int. J. Wildl. Fire https://doi.org/10.1071/WF18091 (2019).
Article  Google Scholar 

47.
Koutsias, N. et al. On the relationships between forest fires and weather conditions in Greece from long-term national observations (1894–2010). Int. J. Wildl. Fire 22, 493–507 (2013).
Google Scholar 

48.
Pereira, M. G., Malamud, B. D., Trigo, R. M. & Alves, P. I. The history and characteristics of the 1980–2005 Portuguese rural fire database. Nat. Hazards Earth Syst. Sci. 11, 3343–3358 (2011).
ADS  Google Scholar 

49.
Belhadj-Khedher, C. et al. A revised historical fire regime analysis in Tunisia (1985–2010) from a critical analysis of the national fire database and remote sensing. Forests 9, 20 (2018).
Google Scholar 

50.
Herrera, S., Bedia, J., Gutiérrez, J. M., Fernández, J. & Moreno, J. M. On the projection of future fire danger conditions with various instantaneous/mean-daily data sources. Clim. Change 118, 827–840 (2013).
ADS  Google Scholar 

51.
Jacob, D. et al. EURO-CORDEX: New high-resolution climate change projections for European impact research. Reg. Environ. Change 14, 563–578 (2014).
Google Scholar 

52.
Moss, R. H. et al. The next generation of scenarios for climate change research and assessment. Nature 463, 747 (2010).
ADS  CAS  PubMed  Google Scholar 

53.
Cannon, A. J. Multivariate quantile mapping bias correction: An N-dimensional probability density function transform for climate model simulations of multiple variables. Clim. Dyn. 50, 31–49 (2018).
Google Scholar 

54.
Cannon, A. J., Sobie, S. R. & Murdock, T. Q. Bias correction of GCM precipitation by quantile mapping: How well do methods preserve changes in quantiles and extremes?. J. Clim. 28, 6938–6959 (2015).
ADS  Google Scholar  More

Our method was implemented in the Google Earth Engine (GEE) platform19. We divided it into four steps. Figure 2 summarizes our approach, including the input of the raw data (land-use and land-cover from 1985 to 2018 and the water surface), and the output data (from 1986 to 2018), which included maps of the annual secondary forest increment (Product 1), annual secondary forest extent (Product 2), annual secondary forest loss (Product 3; from 1987 to 2018), and annual secondary forest age maps (Product 4).
Fig. 2

Workflow of the proposed method.

Full size image

Input data
We used the land-use and land-cover data from the Brazilian Annual Land-Use and Land-Cover Mapping Project (MapBiomas Collection 4.1; https://mapbiomas.org/en/colecoes-mapbiomas-1)1 as input data. This dataset was obtained through the classification of images from the Landsat satellite series (30-m spatial resolution) using a theoretical algorithm implemented in the GEE platform19. Details about the processing of the dataset can be found in the Algorithm Theoretical Basis Document20. More detail about the land-use and land-cover classes can be found in the MapBiomas website (https://mapbiomas.org/en/codigos-de-legenda?cama_set_language=en).
Moreover, we used the maximum water surface extent data (from 1984 to 2018) developed by Pekel et al.21 (https://global-surface-water.appspot.com) to avoid the inclusion of false detection within wetland areas in our products. This dataset contains a map of the spatial distribution of the water surface cover from 1984 to 2018, globally21. These data were obtained from 3,865,618 Landsat 5, 7, and 8 scenes acquired between 16 March 1984 and 31 December 2018. Each pixel was individually classified into water or non-water cover using an expert system21 implemented in the GEE platform19.
Step 1 – Reclassifying MapBiomas data
All MapBiomas land-use and land-cover maps from 1985 to 2018 (34 maps) were reclassified into binary maps. We assigned the value “1” for all pixels in the Forest formation class of the MapBiomas product (Legend ID: 3) and “0” for the other land-use and land-cover classes. In our reclassified maps, pixels with value of “1” were, then, associated to the class “Forest”, which includes only forests classified as old-growth and secondary (before 1985). Mangrove and forest plantation classes were excluded from our secondary forest map.
Step 2 – Mapping the Annual Increment of Secondary Forests
We mapped the annual increment of secondary forests using the forest maps produced in Step 1. This process was carried out pixel-by-pixel, where every pixel classified as Forest (value 1) in the analysed year (yi; between 1986 to 2018) and classified as non-forest (value 0) in the previous year (yi-1; i = 1985, 1986… 2017) was mapped as secondary forest. As forest cover maps before 1985 were not available in the MapBiomas product, maps of secondary forest increment start in 1986, when it was possible to detect the first transition (1985 to 1986). Thus, 33 binary maps were obtained, where the secondary forest increments (non-forest to forest) have a value of 1 and the other transitions a value of 0 (forest to forest, non-forest to non-forest, and forest to non-forest). Here, we only considered secondary forest growth in pixels that had previously an anthropic cover (forest plantation, pasture, agriculture, mosaic of agriculture and pasture, urban infrastructure, and mining) and did not overlap wetland areas.
Step 3 – Mapping the Annual Extent of Secondary Forests
We generated 33 maps of the annual extent of secondary forests. To produce the map of secondary forest extent in 1987, we summed the map of the total secondary forest extent in 1986, which is the same map as the secondary forest increment in 1986 from step 2, with the 1987 increment map, resulting in a map containing all secondary forest pixels from 1986 and 1987. Knowing that the sequential sum of these maps results in pixels with values higher than 1, to create annual binary maps of secondary forest extent, we reclassified the map produced for each year by assigning the value 1 to pixels with values between 2 and 33 (secondary forest extent) and pixels with a value 0 were kept unchanged. Finally, to remove all secondary forest pixels that were deforested in 1987, keeping in the map only pixels with the extent of stand secondary forests, we multiplied the resulting map by the annual forest cover map of 1987, produced in step 1 (Fig. 3). This procedure was applied year-by-year from 1986 to 2018 to produce the maps of annual secondary forest extent. The removal of deforested pixels provides a product depicting the extent of secondary forest deforested in each specific year and they were also included as complimentary maps (from 1987 to 2018) in our dataset.
Fig. 3

Conceptual model of the approach used to calculate the age of secondary forests throughout the Brazilian territory.

Full size image

Step 4 – Calculating the Age of Secondary Forest
Finally, we calculated the age of the secondary forests (Fig. 3). First, we summed the 1986 map of annual secondary forest extent (from Step 3) with the 1987 map to obtain the age of secondary forests in 1987 (Fig. 4). We continued this summation year-by-year until the secondary forest age map of 2018 was obtained (Fig. 4). The values of each pixel in 2018 correspond to the age of the secondary forest. To ensure the elimination of deforested secondary forests from each age map, we executed a similar procedure as described in step 3 by removing all forest pixels overlaying non-forest areas (Fig. 4). As our analyses started in 1986, it was not possible to identify secondary forests before this year. The 1986 age map, therefore, only shows one-year old secondary forests, and the 2018 map shows ages of secondary forest varying between 1 and 33 years (Fig. 4). If a secondary forest pixel with any age is cleared in a given year, it is then removed and a value of zero is attributed to the pixel. The age of this pixel, subsequently, will only be computed again if the algorithm detects a new non-forest to forest transition in the forest cover map (Step 1), which depends on the MapBiomas project classification method.
Fig. 4

(a) Scatter-plot for the relationship between the proportion of the secondary forest within the 10 by 10 km cells in the two datasets. The dashed blue line is the 1:1 line; the red line is the average regression from the bootstrap approach with 10,000 interactions; the dashed red lines are regressions using the standard deviation values of the equation parameters. All p-values from the 10,000 bootstrap interactions were lower than 0.001. (b) Jitter-plot for the proportion of the secondary forest within the 10 by 10 km cells. The red dot is the mean, and the red vertical line the standard deviation.

Full size image More

Affiliations

The Intergovernmental Authority on Development Climate Prediction and Applications Centre (ICPAC), Nairobi, Kenya
Abubakr A. M. Salih, Marta Baraibar, Kenneth Kemucie Mwangi & Guleid Artan

Authors
Abubakr A. M. Salih

Marta Baraibar

Kenneth Kemucie Mwangi

Guleid Artan

Corresponding author
Correspondence to Abubakr A. M. Salih. More

1.
Byrne, M. Impact of ocean warming and ocean acidification on marine invertebrate life history stages: vulnerabilities and potential for persistence in a changing ocean. Oceanogr. Mar. Biol. An Annu. Rev. 49, 1–42 (2011).
Google Scholar 
2.
Pineda, M. C. et al. Tough adults, frail babies: an analysis of stress sensitivity across early life-history stages of widely introduced marine invertebrates. PLoS ONE 7, e46672 (2012).
ADS  CAS  PubMed  PubMed Central  Google Scholar 

3.
Gaines, S. & Roughgarden, J. Larval settlement rate: a leading determinant of structure in an ecological community of the marine intertidal zone. Proc. Natl. Acad. Sci. 82, 3707–3711 (1985).
ADS  CAS  PubMed  Google Scholar 

4.
Pecorino, D., Lamare, M. D., Barker, M. F. & Byrne, M. How does embryonic and larval thermal tolerance contribute to the distribution of the sea urchin Centrostephanus rodgersii (Diadematidae) in New Zealand?. J. Exp. Mar. Bio. Ecol. 445, 120–128 (2013).
Google Scholar 

5.
Pörtner, H. O. & Farrell, A. P. Physiology and climate change. Science 322, 690–692 (2008).
PubMed  Google Scholar 

6.
O’Connor, M. I. et al. Temperature control of larval dispersal and the implications for marine ecology, evolution, and conservation. Proc. Natl. Acad. Sci. U.S.A. 104, 1266–1271 (2007).
ADS  PubMed  PubMed Central  Google Scholar 

7.
Sunday, J. M., Bates, A. E. & Dulvy, N. K. Thermal tolerance and the global redistribution of animals. Nat. Clim. Chang. 2, 686–690 (2012).
ADS  Google Scholar 

8.
Doney, S. C., Fabry, V. J., Feely, R. A. & Kleypas, J. A. Ocean acidification: the other CO2 problem. Ann. Rev. Mar. Sci. 1, 169–192 (2009).
PubMed  Google Scholar 

9.
Fabry, V. J., Seibel, B. A., Feely, R. A., Fabry, J. C. O. & Fabry, V. J. Impacts of ocean acidification on marine fauna and ecosystem processes. ICE J. Mar. Sci. 65, 414–432 (2008).
CAS  Google Scholar 

10.
Kroeker, K. J., Kordas, R. L., Crim, R. N. & Singh, G. G. Meta-analysis reveals negative yet variable effects of ocean acidification on marine organisms. Ecol. Lett. 13, 1419–1434 (2010).
PubMed  Google Scholar 

11.
Beadle, B. Y. L. C. The effect of salinity changes on the water content and respiration of marine invertebrates. J. Exp. Biol. 8, 211–227 (1931).
Google Scholar 

12.
Cheng, B. S., Chang, A. L., Deck, A. & Ferner, M. C. Atmospheric rivers and the mass mortality of wild oysters: Insight into an extreme future?. Proc. R. Soc. B Biol. Sci. 283, 20161462 (2016).
Google Scholar 

13.
Przeslawski, R., Byrne, M. & Mellin, C. A review and meta-analysis of the effects of multiple abiotic stressors on marine embryos and larvae. Global Change Biol. 21, 2122–2140 (2015).
ADS  Google Scholar 

14.
Byrne, M. & Przeslawski, R. Multistressor impacts of warming and acidification of the ocean on marine invertebrates’ life histories. Integ. Comp. Biol. 53, 582–596 (2013).
CAS  Google Scholar 

15.
Kroeker, K. J. et al. Impacts of ocean acidification on marine organisms: Quantifying sensitivities and interaction with warming. Global Change Biol. 19, 1884–1896 (2013).
ADS  Google Scholar 

16.
Bindoff, N. L., et al. Chapter 5: Changing ocean, marine ecosystems, and dependent communities. Intergovernmental panel of climate change. in IPCC Special Report on the Ocean and Cryosphere in a Changing Climate 447–587 (2019).

17.
Feely, R. A. et al. Present and future changes in seawater chemistry due to ocean acidification. Geophys. Monogr. Ser. 183, 175–188 (2009).
CAS  Google Scholar 

18.
Rhein, M. et al. Observations: ocean. In Climate Change 2013 the Physical Science Basis: Working Group I Contribution to the Fifth Assessment Report of the Intergovernmental Panel on Climate Change (eds Stocker, T. F. et al.) (Cambridge University Press, Cambridge, 2013). https://doi.org/10.1017/CBO9781107415324.010.
Google Scholar 

19.
Narita, D., Rehdanz, K. & Tol, R. S. J. Economic costs of ocean acidification: a look into the impacts on global shellfish production. Clim. Change 113, 1049–1063 (2012).
ADS  Google Scholar 

20.
Beck, M. W. et al. Oyster reefs at risk and recommendations for conservation, restoration, and management. Bioscience 61, 107–116 (2011).
Google Scholar 

21.
Blake, B. & Bradbury, A. Washington Department of Fish and Wildlife Plan for Rebuilding Olympia Oyster (Ostrea lurida ) Populations in Puget Sound with a Historical and Contemporary Overview. (2012).

22.
Hettinger, A. et al. The influence of food supply on the response of Olympia oyster larvae to ocean acidification. Biogeosciences 10, 6629–6638 (2013).
ADS  Google Scholar 

23.
Hettinger, A. et al. Larval carry-over effects from ocean acidification persist in the natural environment. Global. Change Biol. 19, 3317–3326 (2013).
Google Scholar 

24.
Li, J. et al. The potential of ocean acidification on suppressing larval development in the Pacific oyster Crassostrea gigas and blood cockle Arcain flata Reeve*. Chin. J. Oceanol. Limnol. 32, 1307–1313 (2014).
ADS  CAS  Google Scholar 

25.
Talmage, S. C. & Gobler, C. J. Effects of elevated temperature and carbon dioxide on the growth and survival of larvae and juveniles of three species of northwest Atlantic bivalves. PLoS ONE 6, e26941 (2011).
ADS  CAS  PubMed  PubMed Central  Google Scholar 

26.
Waldbusser, G. G. et al. Saturation-state sensitivity of marine bivalve larvae to ocean acidification. Nat. Clim. Change 5, 273–280 (2015).
ADS  CAS  Google Scholar 

27.
Dekshenieks, M. M., Hofmann, E. E., Powell, E. N. & Powell, E. N. Environmental effects on the growth and development of eastern Oyster, Crassostrea virginica ( Gmelin, 1791), Larvae : a modeling study. J. Shellfish Res. 12, 241–254 (1993).
Google Scholar 

28.
Ko, G. W. K. et al. Interactive effects of ocean acidification, elevated temperature, and reduced salinity on early-life stages of the pacific oyster. Environ. Sci. Technol. 48, 10079–10088 (2014).
ADS  CAS  PubMed  Google Scholar 

29.
Shanks, A. L., Grantham, B. A. & Carr, M. H. Propagule dispersal distance and the size and spacing of marine reserves. Ecol. Appl. 13, 159–169 (2003).
Google Scholar 

30.
Shanks, A. L. Pelagic larval duration and dispersal distance revisited. Biol. Bull. 216, 373–385 (2009).
PubMed  Google Scholar 

31.
Pineda, J., Hare, J. A. & Sponaugle, S. Larval transport and dispersal in the coastal ocean and consequences for population connectivity. Oceanography 20, 22–39 (2007).
Google Scholar 

32.
Barros, P., Sobral, P., Range, P., Chícharo, L. & Matias, D. Effects of sea-water acidification on fertilization and larval development of the oyster Crassostrea gigas. J. Exp. Mar. Biol. Ecol. 440, 200–206 (2013).
Google Scholar 

33.
Hettinger, A. et al. Persistent carry-over effects of planktonic exposure to ocean acidification in the Olympia oyster. Ecology 93, 2758–2768 (2012).
PubMed  Google Scholar 

34.
Barton, A. et al. Impacts of coastal acidification on the Pacific Northwest shellfish industry and adaptation strategies implemented in response. Oceanography 28, 146–159 (2015).
Google Scholar 

35.
Wasson, K. et al. Coast-wide recruitment dynamics of Olympia oysters reveal limited synchrony and multiple predictors of failure. Ecology https://doi.org/10.1002/ecy.1602 (2016).
Article  PubMed  Google Scholar 

36.
Cole, V. J. et al. Effects of multiple climate change stressors: ocean acidification interacts with warming, hyposalinity, and low food supply on the larvae of the brooding flat oyster Ostrea angasi. Mar. Biol. 163, 1–17 (2016).
CAS  Google Scholar 

37.
Havenhand, J., Dupont, S. & Quinn, G. P. Chapter 4: Designing ocean acidification experiments to maximise inference. in Guide to Best Practices for Ocean Acidification Research and Data Reporting 67–80 (2010).

38.
Mcintyre, B. A., McPhee-Shaw, E. E., Hatch, M. B. & Arellano, S. M. Location matters: passive and active factors affect the vertical distribution of Olympia oyster (Ostrea lurida) larvae. Estuaries Coasts https://doi.org/10.1007/s12237-020-00771-8 (2020).
Article  Google Scholar 

39.
Davis, H. C. On cultivation of larvae of Ostrea lurida. Anat. Rec. 105, 111 (1949).
Google Scholar 

40.
Loosanoff, V. L. & Davis, H. C. Rearing of bivalve mollusks. Adv. Mar. Biol. 1, 1–136 (1963).
Google Scholar 

41.
Hofmann, E. E., Powell, E. N., Bochenek, E. A. & Klinck, J. M. A modelling study of the influence of environment and food supply on survival of Crassostrea gigas larvae. ICES J. Mar. Sci. 61, 596–616 (2004).
Google Scholar 

42.
Barber, J. S., Dexter, J. E., Grossman, S. K., Greiner, C. M. & Mcardle, J. T. Low temperature brooding of Olympia Oysters (Ostrea lurida) in Northern Puget sound. J. Shellfish Res. 35, 351–357 (2016).
Google Scholar 

43.
Hopkins, A. E. Experimental observations on spawining, larval development, and setting in the Olympia oyster Ostrea lurida. Bull U.S.A. Bur. Fish. 48, 439–503 (1937).
Google Scholar 

44.
Pritchard, C., Shanks, A., Rimler, R., Oates, M. & Rumrill, S. The Olympia Oyster Ostrea lurida : recent advances in natural history, ecology, and restoration. J. Shellfish Res. 34, 259–271 (2015).
Google Scholar 

45.
Bible, J. M. et al. Timing of stressors alters interactive effects on a coastal foundation species. Ecology 98, 2468–2478 (2017).
PubMed  Google Scholar 

46.
Waldbusser, G. G. et al. Slow shell building, a possible trait for resistance to the effects of acute ocean acidification. Limnol. Oceanogr. 61, 1969–1983 (2016).
ADS  CAS  Google Scholar 

47.
Lucey, N. M. et al. To brood or not to brood: are marine invertebrates that protect their offspring more resilient to ocean acidification?. Sci. Rep. 5, 1–7 (2015).
Google Scholar 

48.
Barton, A., Hales, B., Waldbusser, G. G., Langdon, C. & Feelyd, R. A. The Pacific oyster, Crassostrea gigas, shows negative correlation to naturally elevated carbon dioxide levels: Implications for near-term ocean acidification effects. Limnol. Oceanogr. 57, 698–710 (2012).
ADS  CAS  Google Scholar 

49.
Miller, A. W., Reynolds, A. C., Sobrino, C. & Riedel, G. F. Shellfish face uncertain future in high CO2 world: influence of acidification on oyster larvae calcification and growth in estuaries. PLoS ONE 4, e5661 (2009).
ADS  PubMed  PubMed Central  Google Scholar 

50.
Khangaonkar, T. et al. Analysis of hypoxia and sensitivity to nutrient pollution in Salish Sea. J. Geophys. Res. Ocean. 123, 4735–4761 (2018).
ADS  CAS  Google Scholar 

51.
Spencer, L. H. et al. Carryover effects of temperature and pCO2 across multiple Olympia oyster populations. Ecol. Appl. 30, e02060 (2020).
PubMed  Google Scholar 

52.
Scheltema, R. S. Larval dispersal as a means of genetic exchange between geographically separated populations of shallow-water benthic marine gastropods. Biol. Bull. 140, 284–322 (1971).
Google Scholar 

53.
Pechenik, J. A. On the advantages and disadvantages of larval stages in benthic marine invertebrate life cycles. Mar. Ecol. Prog. Ser. 177, 269–297 (1999).
ADS  Google Scholar 

54.
Stick, D. A. Identification of Optimal Broodstock for Pacific Northwest Oysters (Oregon State University, Oregon, 2011).
Google Scholar 

55.
Silliman, K. Population structure, genetic connectivity, and adaptation in the Olympia oyster (Ostrea lurida) along the west coast of North America. Evol. Appl. 12, 923–939 (2019).
CAS  PubMed  PubMed Central  Google Scholar 

56.
Hatch, M. B. A., College, N. I. & Wyllie-echeverria, S. Historic distribution of Ostrea lurida (Olympia oyster) in the San Juan Archipelago. Washington State. 1, 38–45 (2008).
Google Scholar 

57.
Polson, M. P. & Zacherl, D. C. Geographic distribution and intertidal population status for the Olympia Oyster, Ostrea lurida Carpenter 1864, from Alaska to Baja. J. Shellfish Res. 28, 69–77 (2009).
Google Scholar 

58.
Dekshenieks, M. M., Hofmann, E. E., Klinck, J. M. & Powell, E. N. Quantifying the effects of environmental change on an Oyster population: a modeling study. Estuaries 23, 593–610 (2000).
Google Scholar 

59.
Powell, E. N., Klinck, J. M., Hofmann, E. E. & Ray, S. M. Modeling oyster populations. IV: Rates of mortality, population crashes, and management. Fish. Bull. 92, 347–373 (1994).
Google Scholar 

60.
Chan, F. et al. Persistent spatial structuring of coastal ocean acidification in the California Current System. Sci. Rep. 7, 1–7 (2017).
Google Scholar 

61.
Long, W. & Khangaonkar, T. Approach for Simulating Acidification and the Carbon Cycle in the Salish Sea to Distinguish Regional Source Impacts. Washington Department of Ecology (2014).

62.
Love, B. A., Olson, M. B. & Wuori, T. Technical note: a minimally invasive experimental system for pCO2 manipulation in plankton cultures using passive gas exchange (atmospheric carbon control simulator). Biogeosciences 14, 2675–2684 (2017).
ADS  CAS  Google Scholar 

63.
Strathmann, M. F. Reproduction and Development of Marine Invertebrates of the Northern Pacific Coast (University of Washington Press, Seattle, 1987).
Google Scholar 

64.
Ko, G. W. K. et al. Larval and post-larval stages of pacific Oyster (Crassostrea gigas) are resistant to elevated CO2. PLoS ONE 8, 1–12 (2013).
Google Scholar 

65.
Buckham, S. Ocean acidification affects larval swimming in Ostrea lurida but not Crassostrea gigas. WWU Graduate School Collection. https://cedar.wwu.edu/wwuet/451 (2015).

66.
Dickson, A., Sabine, C. & Christian, J. (eds). Guide to Best Practices for Ocean CO2 Measurements. In PICES Special Publication 3 191 (2007).

67.
Pelletier, G., Lewis, E. & Wallace, D. co2.sys2.1.xls, a Calculator for the CO2System in Seawater for Microsoft Excel/VBA, Washington State Department of Ecology, Olympia, WA, Brookhaven National Laboratory, Upton, NY. (2012).

68.
Millero, F. J., Graham, T. B., Huang, F., Bustos-Serrano, H. & Pierrot, D. Dissociation constants of carbonic acid in seawater as a function of salinity and temperature. Mar. Chem. 100, 80–94 (2006).
CAS  Google Scholar 

69.
Waldbusser, G. G. et al. Ocean acidification has multiple modes of action on bivalve larvae. PLoS ONE 10, e0128376 (2015).
PubMed  PubMed Central  Google Scholar 

70.
Gazeau, F. et al. Impacts of ocean acidification on marine shelled molluscs. Mar. Biol. 160, 2207–2245 (2013).
CAS  Google Scholar 

71.
Khangaonkar, T., Nugraha, A., Xu, W. & Balaguru, K. Salish Sea response to global climate change, sea level rise, and future nutrient loads. J. Geophys. Res. Ocean. https://doi.org/10.1029/2018JC014670 (2019).
Article  Google Scholar 

72.
Loosanoff, V. L., Davis, H. C. & Chalney, P. E. Dimensions and shapes of larvae of some marine bivalve mollusks. Malacologia 4, 351–435 (1966).
Google Scholar 

73.
Brink, L. A. Molluscs: Bivalvia. Identification Guide to Larval Marine Invertebrates ofthe Pacific Northwest 129–149 (2001).

74.
Hori, J. On the development of the Olympia oyster, Ostrea lurida carpenter, transplanted from United States to Japan. Bull. Jpn. Soc. Sci. Fish 1, 269–276 (1933).
Google Scholar  More

1.
Olesen, J. M., Bascompte, J., Elberling, H. & Jordano, P. Temporal dynamics in a pollination network. Ecology 89, 1573–1582 (2008).
PubMed  Google Scholar 
2.
Petchey, O. L., Brose, U. & Rall, B. C. Predicting the effects of temperature on food web connectance. Philos. Trans. R. Soc. B 365, 2081–2091 (2010).
Google Scholar 

3.
Menke, S., Böhning-Gaese, K. & Schleuning, M. Plant–frugivore networks are less specialized and more robust at forest–farmland edges than in the interior of a tropical forest. Oikos 121, 1553–1566 (2012).
Google Scholar 

4.
Aizen, M. A., Morales, C. L. & Morales, J. M. Invasive mutualists erode native pollination webs. PLoS Biol. 6, e31 (2008).
PubMed  PubMed Central  Google Scholar 

5.
Aizen, M. A., Sabatino, M. & Tylianakis, J. M. Specialization and rarity predict nonrandom loss of interactions from mutualist networks. Science 335, 1486–1489 (2012).
ADS  CAS  PubMed  Google Scholar 

6.
Vázquez, D. P. et al. Species abundance and asymmetric interaction strength in ecological networks. Oikos 116, 1120–1127 (2007).
Google Scholar 

7.
Poisot, T., Stouffer, D. B. & Gravel, D. Beyond species: why ecological interaction networks vary through space and time. Oikos 124, 243–251 (2015).
Google Scholar 

8.
Holt, R. D. & Kotler, B. P. Short-term apparent competition. Am. Nat. 130, 412–430 (1987).
Google Scholar 

9.
May, R. M. Will a large complex system be stable? Nature 238, 413 (1972).
ADS  CAS  PubMed  Google Scholar 

10.
Thébault, E. & Fontaine, C. Stability of ecological communities and the architecture of mutualistic and trophic networks. Science 329, 853–856 (2010).
ADS  PubMed  Google Scholar 

11.
de Ruiter, P. C., Wolters, V., Moore, J. C. & Winemiller, K. O. Food web ecology: playing jenga and beyond. Science 309, 68–71 (2005).
PubMed  Google Scholar 

12.
Ings, T. C. et al. Ecological networks–beyond food webs. J. Anim. Ecol. 78, 253–269 (2009).
PubMed  Google Scholar 

13.
Simanonok, M. P. & Burkle, L. A. Partitioning interaction turnover among alpine pollination networks: spatial, temporal, and environmental patterns. Ecosphere 5, 1–17 (2014).
Google Scholar 

14.
CaraDonna, P. J. et al. Interaction rewiring and the rapid turnover of plant–pollinator networks. Ecol. Lett. 20, 385–394 (2017).
PubMed  Google Scholar 

15.
Petanidou, T., Kallimanis, A. S., Tzanopoulos, J., Sgardelis, S. P. & Pantis, J. D. Long-term observation of a pollination network: fluctuation in species and interactions, relative invariance of network structure and implications for estimates of specialization. Ecol. Lett. 11, 564–575 (2008).
PubMed  Google Scholar 

16.
Kaiser-Bunbury, C. N., Memmott, J. & Müller, C. B. Community structure of pollination webs of mauritian heathland habitats. Perspect. Plant Ecol. Evol. Sys. 11, 241–254 (2009).
Google Scholar 

17.
MacLeod, M., Genung, M. A., Ascher, J. S. & Winfree, R. Measuring partner choice in plant–pollinator networks: using null models to separate rewiring and fidelity from chance. Ecology 97, 2925–2931 (2016).
PubMed  Google Scholar 

18.
Alarcón, R., Waser, N. M. & Ollerton, J. Year-to-year variation in the topology of a plant–pollinator interaction network. Oikos 117, 1796–1807 (2008).
Google Scholar 

19.
Ponisio, L. C., Gaiarsa, M. P. & Kremen, C. Opportunistic attachment assembles plant–pollinator networks. Ecol. Lett. 20, 1261–1272 (2017).
PubMed  Google Scholar 

20.
Burkle, L. A., Marlin, J. C. & Knight, T. M. Plant-pollinator interactions over 120 years: loss of species, co-occurrence, and function. Science 339, 1611–1615 (2013).
ADS  CAS  PubMed  Google Scholar 

21.
Cirtwill, A. R., Roslin, T., Rasmussen, C., Olesen, J. M. & Stouffer, D. B. Between-year changes in community composition shape species roles in an arctic plant–pollinator network. Oikos 127, 1163–1176 (2018).
Google Scholar 

22.
Bascompte, J. & Stouffer, D. B. The assembly and disassembly of ecological networks. Philos. Trans. R. Soc. B 364, 1781–1787 (2009).
Google Scholar 

23.
Jordano, P., Bascompte, J. & Olesen, J. M. Invariant properties in coevolutionary networks of plant–animal interactions. Ecol. Lett. 6, 69–81 (2003).
Google Scholar 

24.
Díaz-Castelazo, C. et al. Changes of a mutualistic network over time: reanalysis over a 10-year period. Ecology 91, 793–801 (2010).
PubMed  Google Scholar 

25.
Tylianakis, J. M., Martínez-García, L. B., Richardson, S. J., Peltzer, D. A. & Dickie, I. A. Symmetric assembly and disassembly processes in an ecological network. Ecol. Lett. 21, 896–904 (2018).
PubMed  Google Scholar 

26.
Gravel, D., Massol, F., Canard, E., Mouillot, D. & Mouquet, N. Trophic theory of island biogeography. Ecol. Lett. 14, 1010–1016 (2011).
PubMed  Google Scholar 

27.
Dáttilo, W., Guimarães, P. R. & Izzo, T. J. Spatial structure of ant–plant mutualistic networks. Oikos 122, 1643–1648 (2013).
Google Scholar 

28.
Poisot, T., Canard, E., Mouillot, D., Mouquet, N. & Gravel, D. The dissimilarity of species interaction networks. Ecol. Lett. 15, 1353–1361 (2012).
PubMed  Google Scholar 

29.
Bramon Mora, B., Gravel, D., Gilarranz, L. J., Poisot, T. & Stouffer, D. B. Identifying a common backbone of interactions underlying food webs from different ecosystems. Nat. Commun. 9, 2603 (2018).
ADS  PubMed  PubMed Central  Google Scholar 

30.
Stouffer, D. B., Sales-Pardo, M., Sirer, M. I. & Bascompte, J. Evolutionary conservation of species roles in food webs. Science 335, 1489–1492 (2012).
ADS  MathSciNet  CAS  PubMed  MATH  Google Scholar 

31.
Baker, N. J., Kaartinen, R., Roslin, T. & Stouffer, D. B. Species roles in food webs show fidelity across a highly variable oak forest. Ecography 38, 130–139 (2015).
Google Scholar 

32.
CaraDonna, P. J. & Waser, N. M. Temporal flexibility in the structure of plant–pollinator interaction networks. Oikos https://doi.org/10.1111/oik.07526 (2020).

33.
Dunne, J. A., Williams, R. J. & Martinez, N. D. Food-web structure and network theory: the role of connectance and size. Proc. Natl Acad. Sci. USA 99, 12917–12922 (2002).
ADS  CAS  PubMed  Google Scholar 

34.
Bascompte, J., Jordano, P., Melián, C. J. & Olesen, J. M. The nested assembly of plant–animal mutualistic networks. Proc. Natl Acad. Sci. USA 100, 9383–9387 (2003).
ADS  CAS  PubMed  Google Scholar 

35.
Chacoff, N. P., Resasco, J. & Vázquez, D. P. Interaction frequency, network position, and the temporal persistence of interactions in a plant–pollinator network. Ecology 99, 21–28 (2018).
PubMed  Google Scholar 

36.
Bastolla, U. et al. The architecture of mutualistic networks minimizes competition and increases biodiversity. Nature 458, 1018 (2009).
ADS  CAS  PubMed  Google Scholar 

37.
Thompson, R. M. et al. Food webs: reconciling the structure and function of biodiversity. Trends Ecol. Evol. 27, 689–697 (2012).
PubMed  Google Scholar 

38.
Tylianakis, J. M., Didham, R. K., Bascompte, J. & Wardle, D. A. Global change and species interactions in terrestrial ecosystems. Ecol. Lett. 11, 1351–1363 (2008).
PubMed  Google Scholar 

39.
Garibaldi, L. A. et al. Stability of pollination services decreases with isolation from natural areas despite honey bee visits. Ecol. Lett. 14, 1062–1072 (2011).
PubMed  Google Scholar 

40.
Goldwasser, L. & Roughgarden, J. Sampling effects and the estimation of food-web properties. Ecology 78, 41–54 (1997).
Google Scholar 

41.
Westphal, C., Steffan-Dewenter, I. & Tscharntke, T. Mass flowering crops enhance pollinator densities at a landscape scale. Ecol. Lett. 6, 961–965 (2003).
Google Scholar 

42.
Miele, V., Ramos-Jiliberto, R. & Vázquez, D. P. Core–periphery dynamics in a plant–pollinator network. Preprint at https://doi.org/10.1101/543637 (2019).

43.
Hackett, T. D. et al. Reshaping our understanding of species’ roles in landscape-scale networks. Ecol. Lett. 22, 1367–1377 (2019).
PubMed  Google Scholar 

44.
Schwarz, B. et al. Temporal scale-dependence of plant–pollinator networks. Oikos https://doi.org/10.1111/oik.07303 (2020).

45.
Bascompte, J. & Melián, C. J. Simple trophic modules for complex food webs. Ecology 86, 2868–2873 (2005).
Google Scholar 

46.
Kondoh, M. Building trophic modules into a persistent food web. Proc. Natl Acad. Sci. USA 105, 16631–16635 (2008).
ADS  CAS  PubMed  Google Scholar 

47.
Vázquez, D. P., Blüthgen, N., Cagnolo, L. & Chacoff, N. P. Uniting pattern and process in plant–animal mutualistic networks: a review. Ann. Bot. 103, 1445–1457 (2009).
PubMed  PubMed Central  Google Scholar 

48.
Cagnolo, L., Salvo, A. & Valladares, G. Network topology: patterns and mechanisms in plant-herbivore and host-parasitoid food webs. J. Anim. Ecol. 80, 342–351 (2011).
PubMed  Google Scholar 

49.
Aizen, M. A. et al. The phylogenetic structure of plant–pollinator networks increases with habitat size and isolation. Ecol. Lett. 19, 29–36 (2016).
PubMed  Google Scholar 

50.
Junker, R. R., Höcherl, N. & Blüthgen, N. Responses to olfactory signals reflect network structure of flower-visitor interactions. J. Anim. Ecol. 79, 818–823 (2010).
PubMed  Google Scholar 

51.
Coux, C., Rader, R., Bartomeus, I. & Tylianakis, J. M. Linking species functional roles to their network roles. Ecol. Lett. 19, 762–770 (2016).
PubMed  Google Scholar 

52.
Bartomeus, I. et al. A common framework for identifying linkage rules across different types of interactions. Funct. Ecol. 30, 1894–1903 (2016).
Google Scholar 

53.
Weinstein, B. G. & Graham, C. H. Persistent bill and corolla matching despite shifting temporal resources in tropical hummingbird-plant interactions. Ecol. Lett. 20, 326–335 (2017).
PubMed  Google Scholar 

54.
Weinstein, B. G. & Graham, C. H. On comparing traits and abundance for predicting species interactions with imperfect detection. Food Webs 11, 17–25 (2017).
Google Scholar 

55.
Eklöf, A. et al. The dimensionality of ecological networks. Ecol. Lett. 16, 577–583 (2013).
PubMed  Google Scholar 

56.
Olito, C. & Fox, J. W. Species traits and abundances predict metrics of plant–pollinator network structure, but not pairwise interactions. Oikos 124, 428–436 (2015).
Google Scholar 

57.
Hart, D. R., Stone, L. & Berman, T. Seasonal dynamics of the lake kinneret food web: the importance of the microbial loop. Limnol. Oceanogr. 45, 350–361 (2000).
ADS  CAS  Google Scholar 

58.
Pilosof, S., Fortuna, M. A., Vinarski, M. V., Korallo-Vinarskaya, N. P. & Krasnov, B. R. Temporal dynamics of direct reciprocal and indirect effects in a host–parasite network. J. Anim. Ecol. 82, 987–996 (2013).
PubMed  Google Scholar 

59.
Holme, P. & Saramäki, J. Temporal networks. Phys. Rep. 519, 97–125 (2012).
ADS  Google Scholar 

60.
Tylianakis, J. M. & Morris, R. J. Ecological networks across environmental gradients. Annu. Rev. Ecol. Evol. Syst. 48, 25–48 (2017).
Google Scholar 

61.
CaraDonna, P. J. Temporal variation in plant-pollinator interactions, Rocky Mountain Biological Laboratory, CO, USA, 2013 – 2015 ver 1. Environmental Data Initiative, https://doi.org/10.6073/pasta/27dc02fe1655e3896f20326fed5cb95f (2020).

62.
Milo, R. et al. Network motifs: simple building blocks of complex networks. Science 298, 824–827 (2002).
ADS  CAS  PubMed  Google Scholar 

63.
Bramon Mora, B., Cirtwill, A. R. & Stouffer, D. B. pymfinder: a tool for the motif analysis of binary and quantitative complex networks. Preprint at https://doi.org/10.1101/364703 (2018).

64.
Pons, P. & Latapy, M. Computing communities in large networks using random walks. J. Graph Algorithms Appl. 10, 191–218 (2006).
MathSciNet  MATH  Google Scholar 

65.
Danon, L., Diaz-Guilera, A., Duch, J. & Arenas, A. Comparing community structure identification. J. Stat. Mech.: Theory E 2005, P09008 (2005).
MATH  Google Scholar 

66.
Koster, J. & McElreath, R. Multinomial analysis of behavior: statistical methods. Behav. Ecol. Sociobiol. 71, 138 (2017).
PubMed  PubMed Central  Google Scholar 

67.
McElreath, R. Statistical Rethinking: A Bayesian Course with Examples in R and Stan (Chapman and Hall/CRC, London, 2018).

68.
Team, S. D. et al. RStan: the R interface to Stan (The R Foundation, 2019). More

Microorganisms and medium composition
C. acetobutylicum DSM 792 and E. aerogenes DSM 30053 were used for all experiments. A modified Clostridium-specific medium without yeast extract was used for growth of mono-culture C. acetobutylicum as previously described in detail elsewhere71. The medium was prepared containing (per L): 0.5 g of KH2PO4, 0.5 g of K2HPO4 and 2.2 g of NH4CH3COO and glucose or cellobiose were added at a concentration of 999 C-mmol. The pH was arranged with 1 mol L−1 NaOH to 6.8. Trace elements solution was prepared as stock 100× solution containing (per L): 0.2 g of MgSO4·7 H2O, 0.01 g of MnSO4·7H2O, 0.01 g of FeSO4·7H2O, 0.01 g of NaCl. Vitamin solution was prepared as stock 200× solution containing (per L): 0.9 g of thiamine, 0.002 g of biotin and 0.2 g of 4-aminobenzoic acid. The trace elements solution and the vitamin solution were used for all experiments. Mono-culture of E. aerogenes was grown in a defined Enterobacter-specific medium, as described elsewhere72. The Enterobacter-specific medium was prepared containing (per L): 13.3 g K2HPO4, 4 g (NH4)2HPO4, 8 mg EDTA and trace elements (2.5 mg CoCl2·6H2O, 15 mg MnCl2·4H2O, 1.5 g CuCl2·4H2O; 3 mg H3BO3; 2.5 mg Na2MoO4·2H2O, 13 mg of Zn(CH3COO)2·2H2O). Glucose and cellobiose were prepared as stock solutions. Media, trace element solution, glucose and cellobiose solutions were flushed with sterile N2 to make the solutions anaerobic and sterilized separately at 121 °C for 20 min. Sterile anaerobic solutions of glucose or cellobiose, trace elements solution and filter sterilized vitamin solution were added into the media before the inoculation inside the sterilized biological safety cabinet (BH-EN 2005, Faster Srl, Ferrara, Italy).
Design of experiments
A mutual medium accommodating the nutritional requirements of both organisms was designed by using the DoE approach. The buffer compositions of two species specific media described above were analysed and the optimum concentrations of AC (NH4Cl), SA (Na+ acetate) and PB (KH2PO4/K2HPO4) capacity were investigated. The setting of DoE for concentration effect of AC, SA and PB capacity was based on 29 randomized runs within concentration range from 3–30 mmol L−1 of AC, 3–150 mmol L−1 of KH2PO4 and 10–120 mmol L−1 of SA (Table 1). Each experiment was performed in triplicates (n = 3), except for set E of the DoE experiment (centre points), which were performed in pentaplicate (n = 5). The DoE experiments were performed twice (N = 2). The end of the exponential growth phase of E. aerogenes and C. acetobutylicum was reached at 45 and 51.5 h, respectively. For modelling, these time points were used. The reason for providing an acetate source in the medium was due to the possibility to add an acetate oxidizing microorganism to the co-culture consortium, which was not performed in the context of this study.
Closed batch cultivations
Cultures of E. aerogenes and C. acetobutylicum were grown anaerobically at 0.3 bar in a 100 Vol.-% N2 atmosphere in a closed batch set-up33. Mono-culture and consortium closed batch experiments were conducted with the final volume of 50 mL medium in 120 mL serum bottles (Ochs Glasgerätebau, Langerwehe, Germany). Each serum bottle contained 45 mL Clostridium-specific medium, Enterobacter-specific medium or E-medium, 0.25 mL vitamin solution, 3.0 mL glucose or cellobiose stock solution, 0.5 mL trace elements solution and 1.25 mL inoculum. The serum bottles were sealed with rubber stoppers (20 mm butyl ruber, Chemglass Life Science LLC, Vineland, USA). For consortium experiments, different inoculum ratios were tested and initial cell concentrations were arranged with the ratios of (E. aerogenes : C. acetobutylicum) 1:2, 1:10, 1:100, 1:1000, 1:10,000 and 1:100,000 at a temperature of 37 °C. Pre-culture of E. aerogenes was diluted in DoE E-medium (Table 1) to inoculate the organism at cell densities of aforementioned ratios. The pressure in the headspace of the serum bottles were measured individually using a manometer (digital manometer LEO1-Ei,−1…3 bar, Keller, Germany). After each measurement, the pressure was released completely from the headspace of serum bottle by penetrating the butyl rubber stopper with a sterile needle. The pressure values were added up to reveal total produced pressure (cumulative pressure). Experiments were performed three times (N = 3) and each set was performed in quadruplicates (n = 4).
Cell counting, absorption measurements, DNA extraction and qPCR
A volume of 1 mL of liquid sample was collected by using sterile syringes at regular intervals for monitoring biomass growth by measuring the absorbance (optical density at 600 nm (OD600)) using a spectrophotometer (Beckman Coulter Fullerton, CA, USA). Every sampling operation was done inside the sterilized biological safety cabinet (BH-EN 2005, Faster Srl, Ferrara, Italy).
E. aerogenes and C. acetobutylicum cells were counted using a Nikon Eclipse 50i microscope (Nikon, Amsterdam, Netherlands) at each liquid/biomass sampling point. The samples for cell count were taken from each individual closed batch run using syringes (Soft-Ject, Henke Sass Wolf, Tuttlingen, Germany) and hypodermic needles (Sterican size 14, B. Braun, Melsungen, Germany). Ten microlitres of sample were applied onto a Neubauer improved cell counting chamber (Superior Marienfeld, Lauda-Königshofen, Germany) with a grid depth of 0.1 mm.
DNA for qPCR was extracted from 1 mL culture samples by centrifugation at 4 °C and 13,400 r.p.m. for 30 min. The following steps were applied for DNA extraction; (1) cells were resuspended in pre-warmed (65 °C) 1% sodium dodecyl sulfate (SDS) extraction buffer and (2) transferred to Lysing Matrix E tubes (MP Biomedicals, Santa Ana, CA, USA) containing an equal volume of phenol/chloroform/isoamylalcohol (25:24:1). (3) Cell lysis was performed in a FastPrep-24 (MP Biomedicals, NY, USA) device with speed setting 4 for 30 s and the lysate was centrifuged at 13,400 r.p.m. for 10 min. (4) An equal volume of chloroform/isoamylalcohol (24:1) was added to the supernatant of the lysate, followed by centrifugation at 13,400 r.p.m. for 10 min and collection of the aqueous phase. (5) Nucleic acids were precipitated with double volume of polyethylenglycol (PEG) solution (30% PEG, 1.6 mol L−1 NaCl) and 1 μL glycogen (20 mg mL−1) as carrier, incubated for 2 h at room temperature. (6) Following centrifugation at 13,400 r.p.m. for 1 h, nucleic acid pellets were washed with 1 mL cold 70% ethanol, dried at 30 °C using a SpeedVac centrifuge (Thermo Scientific, Dreieich, Germany), eluted in Tris-EDTA buffer and stored at −20 °C until further analysis. Nucleic acid quantification was performed with NanoDrop ND-1000 spectrophotometer (NanoDrop Technologies, Wilmington, DE, USA). qPCR assays were developed for quantifying E. aerogenes and C. acetobutylicum in consortium. The primer pairs were designed by targeting species specific genes (Supplementary Table 6) to prevent false-positive amplification and sequences of genes were compared for identifying optimal primer using the ClustalW2 multiple sequence alignment programme (http://www.ebi.ac.uk/Tools/clustalw2/). qPCR assays were performed in Eppendorf Mastercycler epgradientS realplex2 (Eppendorf, Hamburg, Germany). The PCR mixture (20 μL) contained 10 μL SYBR Green labelled Luna Universal qPCR Master Mix (M3003L, New England Biolabs), 0.5 μL of forward and 0.5 μL reverse primer, 8 μL sterile DEPC water and 1 μL of DNA template. Negative controls containing sterile diethyl pyrocarbonate (DEPC) water as a replacement for the DNA templates and DNA template of the non-targeted species were included separately in each run. The amplification protocol started with an initial denaturation at 95 °C for 2 min, followed by 45 cycles of denaturation at 95 °C for 30 s, annealing and fluorescence acquisition at 60 °C for 30 s and elongation at 72 °C for 30 s. A melting-curve analysis (from 60 °C to 95 °C at a transition rate of 1 °C every 10 s) was performed to determine the specificity of the amplification. All amplification reactions were performed in triplicates. A standard curve was generated as described elsewhere29. Culture samples of each organism were collected at different time intervals for cell count and genomic DNA extraction cell density of each strain were determined by cell counting under microscope during growth and subsequent gDNA extraction was applied to reflect absolute quantification. Six tenfold dilution standards were prepared and a linear regression analysis was performed between qPCR reads and cell counts and OD600 measurements.
Quantification of gas composition
Gas chromatography (GC) measurements were performed from serum bottles that remained without any manipulation after inoculation until the first time point GC measurement. After every GC measurement, remaining gas was released completely from the serum bottles by penetrating the butyl rubber stopper using a sterile needle. The pressure of serum bottles headspace was determined to examine whether there was any remaining overpressure by using a manometer (digital manometer LEO1-Ei,−1…3 bar, Keller, Germany). The gas compositions were analysed by using a GC (7890 A GC System, Agilent Technologies, Santa Clara, USA) with a 19808 Shin Carbon ST Micropacked Column (Restek GmbH, Bad Homburg, Germany) and provided with a gas injection and control unit (Joint Analytical System GmbH, Moers, Germany) as described before73,74,75. The standard test gas employed in GC comprised the following composition: 0.01 Vol.-% CH4; 0.08 Vol.-% CO2 in N2 (Messer GmbH, Wien, Austria). All chemicals were of highest grade available. H2, CO2, N2, 20 Vol.-% H2 in CO2 and 20 Vol.-% CO2 in N2 were of test gas quality (Air Liquide, Schwechat, Austria).
Quantification of liquid metabolites
Quantification of sugars, volatile fatty acids and alcohols were performed with high-performance liquid chromatography (HPLC) system (Agilent 1100), consisting of a G1310A isocratic pump, a G1313A ALS autosampler, a Transgenomic ICSep ICE-ION-300 column, a G1316A column thermostat set at 45 °C and a G1362A RID refractive index detector, measuring at 45 °C (all modules were from Agilent 1100 (Agilent Technologies, CA, USA). The measurement was performed with 0.005 mol L−1 H2SO4 as solvent, with a flow rate of 0.325 mL min−1 and a pressure of 48–49 bar. The injection volume was 40 µL.
Data analysis
For the quantitative analysis, the maximum specific growth rate (µmax [h−1]) and mean specific growth rate (µmean [h−1]) were calculated as follows: N = N0·eµt with N, cell number [cells ml−1]; N0, initial cell number [cells ml−1]; t, time [h] and e, Euler’s number. According to the delta cell counts in between sample points, µ was assessed. The Y(H2/S) [mol mol−1], HER [mmol L−1 h−1], CER [mmol L−1 h−1] and the specific H2 production rate (qH2) [mmol g−1 h−1]32 were calculated from the intervals between each time point and the gas composition in the headspace of serum bottle was determined using the GC. The elementary composition of the corresponding biomass59 was used for the calculation of the mean molar weight, carbon balance and the DoR balance. Yields of byproducts were determined after HPLC measurement. Values were normalized according to the zero control. Moreover, the Shannon diversity index (H) was calculated to interpret the changes in microbial diversity, accounting for both richness (S), the number of species present and abundance of different species. Relative abundance of two species was evaluated according to the calculated evenness (EH) values76. Global substrate uptake rate, byproduct production rates and the mass balance analyses of the mono-cultures and consortium on glucose and cellobiose were calculated between the first and last time point.
Fluorescence in situ hybridization
For FISH, samples of 2 mL were collected for cell fixation. The samples were centrifuged in micro-centrifuge (5415-R, Eppendorf, Hamburg, Germany) for 10 min at 13,200 r.p.m. and pellets were resuspended in 0.5 mL phosphate-buffered saline (PBS) (10 mmol L−1 of Na2HPO4/NaH2PO, 130 mmol L−1 of NaCl, pH of 7.2–7.4). After repeating this procedure twice, 0.5 mL ice-cold absolute ethanol was added to the 0.5 mL PBS/cell mixture. The ethanol fixed samples were thoroughly mixed and then stored at −20 °C. Poly-l-lysine solution (0.01 % (v/v)) was used for coating the microscope slides (76 × 26 × 1 mm, Marienfeld-Superior, Lauda-Königshofen, Germany) containing ten reaction wells separated by an epoxy layer. After dipping the slide into the solution for 5 min, residual poly-l-lysine from the slides was removed by draining the well, followed by air-drying for several minutes. Cells were immobilized on prepared slides by adding samples (1–10 µL) on each well and air-drying. For cell dehydration, the slides were impregnated with ethanol concentrations of 50% (v/v), 80% (v/v) and 96% (v/v), respectively. The slides were dipped into each solution for 3 min, starting from the lowest concentration.
The EUB338 probe77 was used to target specific 16S rRNA found in almost all organisms belonging to the domain of bacteria78. The GAM42a probe specifically binds to target regions of gammaproteobacterial 23S rRNA79 (Supplementary Table 7). Both probes were diluted with DEPC water to a certain extent depending on the fluorescence label. Cy3-labelled EUB338 was diluted to a probe concentration of 30 ng DNA μL−1, whereas FLUOS-labelled GAM42a was adjusted to a final concentration of 50 ng DNA μL−1. For hybridization of the probe, 20 µL of hybridization buffer (900 mmol L−1 NaCl, 20 mmol L−1 Tris/HCl, 30% formamide (v/v), 0.01% SDS (v/v)) and 2 µL of diluted probe solution were added into each well. The hybridization reaction (46 °C, overnight) was facilitated using an airtight hybridization chamber (50 mL centrifuge tube) to prevent dehydration.
A stringent washing step was performed at 48 °C for 10 min in pre-warmed 50 mL washing buffer (100 mmol L−1 NaCl, 20 mmol L−1 Tris/HCl, 5 mmol L−1 EDTA). Afterwards, the slides were dried up and a mounting medium (Antifade Mounting Medium, Vectashield Vector Laboratories, CA, USA) was added to each well. The slides were sealed with a cover glass and examined under phase-contrast microscope (Nikon Eclipse Ni equipped with Lumen 200 Fluorescence Illumination Systems) using filter sets TRITC (557/576) (maximum excitation/emission in nm) for cy3-labelled EUB338 probe and FITC (490/525) for FLUOS-labelled GAM42a probes by a 100 × 1.45 numerical aperture microscope objective (CFI Plan Apo Lambda DM ×100 Oil; Nikon Corp., Japan).
Statistics and reproducibility
DoE experiments were designed and analysed using Design Expert version 11.1.2.0 (Stat-Ease, Inc. USA). Analysis of variation was performed at α = 0.05. The p-values for each test are indicated in the ‘Results’ section. All closed batch experiments were reproduced three times (N = 3) and each replication contained quadruplicate (n = 4). qPCR and FISH experiments, which applied all of the mentioned replicates, were performed in technical triplicates (n = 3). DoE experiments were conducted twice (N = 2) and each replication contained triplicate experiments for corner points (n = 3), except the set E (centre points), which was performed in biological pentaplicates (n = 5).
Reporting summary
Further information on research design is available in the Nature Research Reporting Summary linked to this article. More

1.
Ollerton, J., Winfree, R. & Tarrant, S. How many flowering plants are pollinated by animals?. Oikos 120, 321–326 (2011).
Article  Google Scholar 
2.
Klein, A. M. et al. Importance of pollinators in changing landscapes for world crops. Proc. R. Soc. B Biol. Sci. 274, 303–313 (2007).
Article  Google Scholar 

3.
Aizen, M., Garibaldi, L. A., Cunningham, S. & Klein, A. M. Long-term global trends in crop yield and production reveal no current pollination shortage but increasing pollinator dependency. Curr. Biol. 18, 1572–1575 (2008).
CAS  Article  Google Scholar 

4.
Aizen, M., Garibaldi, L. A., Cunningham, S. & Klein, A. M. How much does agriculture depend on pollinators? Lessons from long-term trends in crop production. Ann. Bot. 103, 1579–1588 (2009).
Article  Google Scholar 

5.
Kremen, C., Williams, N. M. & Thorp, R. W. Crop pollination from native bees at risk from agricultural intensification. Proc. Natl. Acad. Sci. USA. 99, 16812–16816 (2002).
ADS  CAS  Article  Google Scholar 

6.
Klein, A. M., Steffan-Dewenter, I. & Tscharntke, T. Fruit set of highland coffee increases with the diversity of pollinating bees. Proc. R. Soc. B Biol. Sci. 270, 955–961 (2003).
Article  Google Scholar 

7.
Hoehn, P., Tscharntke, T., Tylianakis, J. M. & Steffan-Dewenter, I. Functional group diversity of bee pollinators increases crop yield. Proc. R. Soc. B Biol. Sci. 275, 2283–2291 (2008).
Article  Google Scholar 

8.
Garibaldi, L. A. et al. Wild pollinators enhance fruit set of crops regardless of honey bee abundance. Science 339, 1608–1611 (2013).
ADS  CAS  Article  Google Scholar 

9.
Potts, S., Imperatriz-Fonseca, V. & Ngo, H. The assessment report of the Intergovernmental Science-Policy Platform on Biodiversity and Ecosystem Services on pollinators, pollination and food production. (Secretariat of the Intergovernmental Science-Policy Platform on Biodiversity and Ecosystem Services, 2016). https://doi.org/10.5281/zenodo.3402856

10.
Leong, J. M. & Thorp, R. W. Colour-coded sampling: the pan trap colour preferences of oligolectic and nonoligolectic bees associated with a vernal pool plant. Ecol. Entomol. 24, 329–335 (1999).
Article  Google Scholar 

11.
Westphal, C. et al. Measuring bee diversity in different European habitats and biogeographical regions. Ecol. Monogr. 78, 653–671 (2008).
Article  Google Scholar 

12.
Wilson, J. S., Griswold, T. & Messinger, O. J. Sampling bee communities (Hymenoptera: Apiformes) in a desert landscape: are pan traps sufficient?. J. Kansas Entomol. Soc. 81, 288–300 (2008).
Article  Google Scholar 

13.
Toler, T. R., Evans, E. W. & Tepedino, V. J. Pan-trapping for bees (Hymenoptera: Apiformes) in Utah’s west desert: The importance of color diversity. Pan-Pac. Entomol. 81, 103–113 (2005).
Google Scholar 

14.
Nielsen, A. et al. Assessing bee species richness in two Mediterranean communities: importance of habitat type and sampling techniques. Ecol. Res. 26, 969–983 (2011).
Article  Google Scholar 

15.
Saunders, M. E. & Luck, G. W. Pan trap catches of pollinator insects vary with habitat. Aust. J. Entomol. 52, 106–113 (2013).
Article  Google Scholar 

16.
Allen-Wardell, G. et al. The potential consequences of pollinator declines on the conservation of biodiversity and stability of food crop yields. Conserv. Biol. 12, 8–17 (1998).
Article  Google Scholar 

17.
Kearns, C., Inouye, D. & Waser, N. Endangered mutualisms: the conservation of plant-pollinator interactions. Annu. Rev. Ecol. Syst. 29, 83–112 (1998).
Article  Google Scholar 

18.
Brunet, J. Pollinator decline: implications for food security and environment. Sci. Glob. https://doi.org/10.33548/scientia371 (2019).
Article  Google Scholar 

19.
Popic, T. J., Davila, Y. C. & Wardle, G. M. Evaluation of common methods for sampling invertebrate pollinator assemblages: net sampling out-perform pan traps. PLoS ONE 8, e66665 (2013).
ADS  CAS  Article  Google Scholar 

20.
Bauer, A. A., Clayton, M. K. & Brunet, J. Floral traits influencing plant attractiveness to three bee species: consequences for plant reproductive success. Am. J. Bot. 104, 1–10 (2017).
Article  Google Scholar 

21.
Brunet, J. & Stewart, C. M. Impact of bee species and plant density on alfalfa pollination and potential for gene flow. Psyche A J. Entomol. 2010, 1–7 (2010).
Article  Google Scholar 

22.
Wang, X. et al. Biodiversity of wild alfalfa pollinators and their temporal foraging characters in Hexi Corridor Northwest China. Entomol. Fenn. 23, 4–12 (2012).
Article  Google Scholar 

23.
Chen, M., Zhao, X. Y. & Zuo, X. A. Pollinator activity and pollination success of Medicago sativa L. in a natural and a managed population. Ecol. Evol. 8, 9007–9016 (2018).
Article  Google Scholar 

24.
Cane, J. H. Pollinating bees (Hymenoptera: Apiformes) of U.S. alfalfa compared for rates of pod and seed set. J. Econ. Entomol. 95, 22–27 (2002).
Article  Google Scholar 

25.
Bohart, G. E. Alfalfa pollinators with special reference to species other than honey bees. In Proceedings of the 10th International Congress of Entomology, Vol. 4, pp. 929–937 (1958).

26.
Brookes, B., Small, E., Lefkovitch, L. P., Damman, H. & Fairey, D. T. Attractiveness of alfalfa (Medicago sativa L.) to wild pollinators in relation to wildflowers. Can. J. Plant Sci. 74, 779–783 (1994).
Article  Google Scholar 

27.
Bohart, G. E. Pollination of alfalfa and red clover. Annu. Rev. Entomol. 2, 355–380 (1957).
Article  Google Scholar 

28.
Chao, A. et al. Rarefaction and extrapolation with Hill numbers: a framework for sampling and estimation in species diversity studies. Ecol. Monogr. 84, 45–67 (2014).
Article  Google Scholar 

29.
Hall, H. G. Color preferences of bees captured in pan traps. J. Kansas Entomol. Soc. 89, 273–276 (2016).
Article  Google Scholar 

30.
Campbell, J. W. & Hanula, J. L. Efficiency of Malaise traps and colored pan traps for collecting flower visiting insects from three forested ecosystems. J. Insect Conserv. 11, 399–408 (2007).
Article  Google Scholar 

31.
Heneberg, P. & Bogusch, P. To enrich or not to enrich? Are there any benefits of using multiple colors of pan traps when sampling aculeate Hymenoptera?. J. Insect Conserv. 18, 1123–1136 (2014).
Article  Google Scholar 

32.
Moreira, E. F. et al. Are pan traps colors complementary to sample community of potential pollinator insects?. J. Insect Conserv. 20, 583–596 (2016).
Article  Google Scholar 

33.
Burd, M. Bateman’s principle and plant reproduction: the role of pollen limitation in fruit and seed set. Bot. Rev. 60, 83–139 (1994).
MathSciNet  Article  Google Scholar 

34.
Herrera, C. M. Pollinator abundance, morphology, and flower visitation rate: analysis of the ‘quantity’ component in a plant-pollinator system. Oecologia 80, 241–248 (1989).
ADS  Article  Google Scholar 

35.
Riday, H., Reisen, P., Raasch, J. A., Santa-Martinez, E. & Brunet, J. Selfing rate in an alfalfa seed production field pollinated with leafcutter bees. Crop Sci. 55, 1087–1095 (2015).
Article  Google Scholar 

36.
McGregor, S. Insect Pollination of Cultivated Crop Plants. (USDA, 1976). https://doi.org/10.1093/besa/23.1.104

37.
Grundel, R., Frohnapple, K. J., Jean, R. P. & Pavlovic, N. B. Effectiveness of bowl trapping and netting for inventory of a bee community. Environ. Entomol. 40, 374–380 (2011).
Article  Google Scholar 

38.
Oksanen, J. et al. Vegan: community ecology package. R package version 2.5–5. https://CRAN.R-project.org/package=vegan (2019).

39.
R Core Team. R: A language and environment for statistical computing. (R Foundation for Statistical Computing, 2019).

40.
Jost, L. Entropy and diversity. Oikos 113, 363–375 (2006).
Article  Google Scholar 

41.
Hsieh, T. C., Ma, K. H. & Chao, A. iNEXT: an R package for rarefaction and extrapolation of species diversity (Hill numbers). Methods Ecol. Evol. 7, 1451–1456 (2016).
Article  Google Scholar 

42.
Signorell, A. & Al, E. DescTools: tools for descriptive statistics. R package version 0.99.28. (2019). More