The large mean body size of mammalian herbivores explains the productivity paradox during the last glacial maximum

Large herbivores are a major agent in ecosystems, influencing vegetation structure and carbon and nutrient flows. During the last glacial period, the steppe-tundra ecosystem prevailed on the unglaciated northern lands, hosting a high diversity and density of megafaunal herbivores. The apparent discrepancy between abundant megafauna and the expected low vegetation productivity under a generally harsher climate with lower CO2 concentration, termed productivity paradox, awaits large-scale quantitative analysis from process-based ecosystem models. Yet most of the current global dynamic vegetation models (DGVMs) lack explicit representation of large herbivores. Here we incorporated a grazing module in the ORCHIDEE-MICT DGVM based on physiological and demographic equations for wild large grazers, taking into account feedbacks of large grazers on vegetation. The model was applied globally for present-day and the last glacial maximum (LGM). The present-day results of potential grazer biomass, combined with an empirical land use map, infer a reduction of wild grazer biomass by 79-93% due to anthropogenic land replacement over natural grasslands. For the LGM, we find that the larger mean body size of mammalian herbivores than today is the crucial clue to explain the productivity paradox, due to a more efficient exploitation of grass production by grazers with a larger-body size.


Supplementary
. Correlation of body mass and dry matter intake rate for mammalian herbivores. The data were compiled in Clauss et al. 7 Blue circles indicate all the 93 herbivorous mammals compiled in ref 7 , while red circles indicate the subset of 46 large herbivores with body mass larger than 10 kg. The allometric regression equation of the subset of large herbivores gives a higher exponent (0.85). Table 1. Herbivore biomass density of protected areas in Africa and Asia, and ecosystems in North America, calculated from the raw data provided by Hatton et al. 6 The data for another 15 ecosystems in North America 8 (Messier, 1994) cited in Hatton et al. 6 were not used in this study, since they only included one herbivore species (moose). Supplementary Table 2. Present-day (PD) and pre-industrial (PI) global and regional total grazer biomass (unit: million tonnes live weight) simulated by ORCHIDEE-MICT. The model simulates natural vegetation without anthropogenic land use, and thus potential values. The gridded fractional covers of the categories "Wild" and "Seminatural" from Antromes version 2 3 for the year 2000 and 1900 were multiplied to the PD and PI potential values, in order to estimate the reduction of wild large grazers due to human land use, with the relative reductions listed in parenthesis.

Supplementary Note 1. Subtraction of tropical rainforest from simulated potential large grazer density
The current model, lacking an explicit representation of browsers and mixed feeders, produced a significant distribution of potential large grazers in densely forested areas, in particular in the tropical rainforests ( Supplementary Fig. 8a). This is because the grazing module is primarily driven by grass productivity, and the grasses in tropical forested areas, although with a very small simulated fractional land cover (ca. 10%), are highly-productive under the favourable climatic conditions and can still support some grazers (ca. 1000 kg km -2 , Supplementary Fig. 8b) in the model. In reality, large herbivores in tropical rainforests are mainly generalist feeders like elephants, feeding on both grass and browse plants 9 ; whereas conventional grazers do not usually live in the rainforest (with a few exceptions like capybaras in South America). To alleviate this model drawback, we subtracted the fractional covers of tropical rainforest, according to the empirically-derived potential natural vegetation map 10 , from the direct model output of potential grazer density ( Supplementary Fig. 8a), and thus derived Fig. 2a and the total regional/global values listed in Supplementary Table 1. We did not do similar subtractions for boreal and temperate forests, because: 1) these forests are not totally closed due to disturbances from extreme weather events like drought, frost and storms, fires, pathogen outbreaks 11,12 , and potentially, large herbivores. The fractional tree cover in these forests seldom exceeds 80%, according to the remote sensing-derived product of vegetation continuous field 13 ; and 2) some large grazer species can indeed live in the northern forests, such as aurochs 14 , the extinction of which has been attributed to human hunting and expulsion 15 . Still, our model may have overestimated grazer density in nontropical dense forests, because most large herbivore species inhabiting forests are mixed feeders and browsers, which are missing in the current model.
For the LGM results, no such subtraction was applied either, considering the very limited distribution of dense forests (e.g. areas with modelled tree cover > 0.8 accounted for only 6% of total land area for the tropics, in contrast to 42% in the present-day simulation). Figure 8. (a) Modelled potential large grazer biomass density for presentday (1960-2009 mean), without post-processing. (b) Relationship between modelled grazer biomass density and tree cover for present-day. Grid cells in tropical and extratropical regions were divided into tree cover bins of 0.04, with the circles representing median values and the error bars indicating 25 th -75 th percentiles for each bin. The size of each circle is proportionate to the number of pixels in each bin.

Supplementary Note 2. Conversion from the modelled plant functional types (PFTs) into mega-biomes
In order to facilitate comparison between modelled LGM vegetation distribution and reconstructions based on pollen and plant macrofossil data from BIOME 6000 version 4.2 (available online: http://www.bridge.bris.ac.uk/resources/Databases/BIOMES_data), we used a method similar to that of Prentice et al. 16 and Kageyama et al. 17 to convert modelled vegetation properties into the 9 "mega-biomes" provided by BIOME 6000. The algorithm is shown in Supplementary Fig. 9: firstly, the global vegetation is divided between cold biomes and the rest, with the threshold GDD5 (annual growing degree days above 5 ºC) equalling to 350 K days; secondly, the separation between forest/savanna, dry grassland/shrubland, and desert (or between forest/dry tundra and tundra) depends on the modelled total foliage projective cover (FPC), which is a function of both fractional cover and leaf area index (LAI) of each PFT ( ) 1 ( ); thirdly, forest and savanna (or forest and dry tundra) are separated by the average height of all existing tree PFTs; finally, within the forest, the dominance of particular tree PFTs is used to separate tropical, temperate and boreal forests. After the conversion, only one biome is assigned to each grid cell. This conversion algorithm makes use of not only modelled PFT fractional covers, but also FPC and average height, which relates to modelled LAI and woody biomass and ultimately relates to vegetation productivity and allocation in the model. Thus it enables the distinction between the relatively more productive dry tundra (corresponding to graminoid and forb tundra, see Table 3 in Harrison and Prentice 18 ) and the low productive tundra (corresponding to cushion forb tundra, erect dwarf shrub tundra, etc.) in Siberia and Alaska. Figure 9. Algorithm to convert the modelled PFT properties into the 9 megabiomes provided by BIOME 6000, adapted from Prentice et al. 16 and Kageyama et al. 17 Supplementary discussion. Impacts of grazing on land carbon cycle during the LGM.

Supplementary
As shown in Fig. 5, total aboveground grass consumption by grazers was 1.7 Pg C yr -1 , or 3.2% of the total body weight in the unit of dry matter intake per day, which is within the empirical range for domestic livestock (2-4%) 19 . Modelled global mean ratio of grazer-tograss biomass was 0.5% (in g C m -2 : g C m -2 ), comparable to the estimate of 0.8% for the ratio of terrestrial herbivore-to-primary producer biomass by Harfoot et al. 20 using the Madingley ecosystem model, knowing that their result included all plant-eating animals, not only the large mammalian grazers considered in this study. Our result is also within the range of an empirical estimate of 1.1±0.7% for the herbivore-to-primary producer biomass ratio in temperate and tropical grasslands 21 . Figure 10. Impacts of grazing on tree cover and grassland productivity at the LGM. (a) Difference of tree fractional cover between LGM-withGrazer and LGM-noGrazer.

Supplementary
(b) Ratio of grass NPP (gC m -2 grassland yr -1 ) between LGM-withGrazer and LGM-noGrazer. In order to separate the effect of changing grass fractional cover, we plotted grass NPP per unit area of grass instead of per unit area of ground ("grass NPP" elsewhere in this study refers to the latter). Supplementary Fig. 10 displays the difference in tree cover and grassland productivity between model results with and without grazers. The difference in tree cover showed a general reduction with grazing ( Supplementary Fig. 10a), mainly due to the tramplinginduced tree mortality (equation (10) in Methods). In the model, grazers also have an indirect positive impact on the trees through the interaction between fire and vegetation, because grazing reduces grass fuel load and thus the frequency of fires 22,23 , which leads to smaller fireinduced tree mortalities. The slight increase of tree fractional cover in few grid cells ( Supplementary Fig. 10a) corresponded to the areas where grazers favour trees by causing a reduction of fire occurrence that exceeds the negative effects of trampling in the model.
The suppression of woody plants by large grazers and the subsequent alteration in landscape structure has been observed in modern exclosure experiments and paleoecological records 24 . A remote-sensing observation in Kruger National Park in South Africa shows an average change of woody cover (vegetation > 1 m tall) from 20% to 12%, or 40% reduction, in areas accessible to herbivores, compared to the areas of long-term herbivore exclusions 25 . Among the herbivores with different body size, elephants have been revealed as the primary agent of treefall in savannas, due to their physical strength and height causing disproportionate mortality of shrubs and maturing trees 24,26 ; while smaller herbivores, especially browsers, increase mortalities of tree saplings and suppress tree regeneration 27 . Our model results show a general reduction of 5-10% absolute tree fractional cover in tropical and sub-tropical regions with the simulated grazer density during the LGM (Supplementary Fig. 10a), close to the observed reduction in ref 25 . In contrast, simulations with the LPJ-GUESS model over Africa with a parameterization of wild grazers 23 showed no significant changes in vegetation distribution with and without grazers. This is partly because of the different parameterizations of trampling effect, as we calibrated the trampling-induced tree mortality (equation (10)) to match the mortality caused by elephants 28 , while in LPJ-GUESS, which did not yet include a herbivore functional type to represent elephants, the trampled area by grazers was used to calculate the decrease of tree saplings 23 , which may underestimate the disproportionate effect of elephants compared to the one of smaller wild herbivores 26 .
Grazing increased grass NPP per unit area of grass in most of the grid cells ( Supplementary  Fig. 10b), with larger increments in areas of higher grazer biomass density, indicating that the positive effects of grazers on grass NPP in the model, i.e. the regrowth of more productive young leaves after defoliation and the higher photosynthesis capacities (a process simply calibrated from herbivore exclusion experiments to approximate the positive effects associated with accelerated nutrient cycling and traits/composition changes -see Methods), outweighed the negative effect of biomass removal by grazing. The above-to belowground ratio of global grass NPP increased from 1.2 in the LGM without grazers to 1.3 with grazers, due to more allocation of assimilated carbon to leaves to compensate for defoliation in our model 29 . Such a preferential allocation to the shoot system compared to roots in response to defoliation is supported by experimental evidence 30 . As for grass standing biomass, total aboveground biomass was slightly smaller with grazers (5.2 Pg C) than without grazers (5.4 Pg C), resulting from intake being not fully compensated by a higher NPP; while belowground biomass increased from 4.2 Pg C without grazers to 6.6 Pg C with grazers due to a higher NPP.