Human-induced nutrient eutrophication is a major threat to grassland biodiversity, because it promotes the dominance of fast-growing plants. Negative impacts of fertilization on plant biodiversity may be offset by grazing by large vertebrate herbivores. However, whether grazers also mitigate impacts of nutrient addition on insects is less well understood. We use a field experiment to test how plant communities and abundances of pollinators and grasshoppers respond to nutrient addition and grazing by different assemblages of large herbivores, i.e., access by all herbivores (including cattle and horses), access by wild herbivores only (wild boar and deer), no access by large herbivores. Plant biomass increased, plant diversity decreased and community composition shifted towards lower forb cover in response to fertilization, but only in the absence of all herbivores. Flower visitation by Hymenoptera (bees and wasps), i.e., the most abundant pollinator group, was reduced by nutrient addition only in the absence of all herbivores and was positively related to flowering plant richness. In contrast, flower visitation by Diptera (e.g., hoverflies) was enhanced by fertilization, but not affected by grazing. Orthoptera (grasshopper) abundance was reduced by grazing and enhanced by nutrient addition, with positive impacts of fertilization tending to be stronger in plots with only wild or no herbivores. The abundance of grasshoppers was positively related to grass biomass. We conclude that vertebrate herbivores can offset impacts of fertilization on both plant and insect communities, making grazing by large mammals an essential tool to protect insects, particularly pollinators. Most responses to nutrient addition were only apparent in plots without any large herbivores, suggesting that wild herbivores alone could already mitigate nutrient impacts. We also show that insects with contrasting feeding guilds may be favoured by fertilized, ungrazed conditions. Therefore, creating a mosaic of patches grazed at different intensities will enhance overall insect biodiversity.

Study system: We set up a field experiment with nutrient addition and grazing treatments in a semi-natural grassland in the Veluwe area in the Netherlands (52.03 °N, 5.80 °E) (e.g., Kardol et al. 2006, Hannula et al. 2017, Veen et al. 2018). The field site is located on a loamy, sandy mineral soil (Kardol et al. 2006). The mean annual temperature was around 10.7°C and the mean annual precipitation was approximately 840 mm (Veen et al., 2018; Royal Netherlands Meteorological Institute (KNMI)). The field has been used for agriculture and was abandoned in 1982. Since then, it has been managed by Natuurmonumenten as nature reserve and has developed into a grassland with a Plantagini-Festucion plant community (Kardol et al. 2006). The grassland is visited by naturally occurring herbivores, including roe deer (Capreolus capreolus), mainly a browser, and red deer (C. elaphus), a mixed feeder, and wild boar (S. scrofa), a facultative herbivore. European hares and rabbits occur in the area, but at very low densities. In addition, domesticated semi-wild herbivores, i.e., Sayaguesa cattle (Bos taurus) and wild horses (Equus ferus), both grazers, have been introduced as proxies for their extinct ancestor’s aurochs and tarpan and roam year-round freely in the area. In the seven years prior to our study, the densities of cattle and horses were both fluctuating around 2.5 animals (~1875 kg and 1050 kg, respectively) per 100 hectares; the densities of red deer and wild boar varied between 5-7.5 animals (i.e., ~110-170 kg) and 2.5-11.5 animals (i.e., ~190-860 kg) per 100 hectares, respectively. In general, densities of large grazers in Europe are relatively low as a result of management, including hunting; however in rewilding areas, such as our study site, densities of herbivores are generally higher, better resembling historical densities (Fløjgaard et al. 2022).

Experimental design: In the fall of 2017, we established a controlled nutrient addition × grazing experiment, which is part of the Nutrient Network (NutNet) (Borer et al. 2014a). In our experiment, we established five independent, randomized blocks. Each block consisted of six plots of 5 × 5 m, which were grouped in three pairs (Fig. S1). One pair of plots was fenced with 2 m high fences to exclude all large, vertebrate herbivores (referred to as “exclosure”); one pair of plots was fenced with barbed wires allowing access to wild boar and deer while excluding domesticated herbivores (referred to as “wild herbivores”); and one pair of plots was left unfenced to allow grazing by all herbivores (referred to as “all herbivores”). Within each pair, one plot has been fertilized with slow-release grains containing nitrogen (N; Cote N (2M) 43% urea,), phosphorous (P; Triplesuperphosphate, 45% P₂O₅) and potassium (K, Potassiumsulfate, K₂SO₄) each year starting in 2018 at the start of the growing season (end of March) according to Borer et al. (2014a). For N, P, and K we apply 10 g of each nutrient per m^-2 yr^-1. In 2018, at the start of the experiment, all fertilized plots also received a mix of micronutrients, equalling 100 g m^-2 and including Fe (15%), S (14%), Mg (1.5%), Mn (2.5%), Cu (1%), Zn (1%), B (0.2%) and Mo (0.05%) (Borer et al. 2014a). All fertilizers were ordered at Telermaat, Boskoop, the Netherlands.

Data collection: We measured aboveground plant biomass and plant community composition at the peak of the growing season in July 2017 (starting conditions before establishing exclosures) and after three years of fertilization and herbivore exclusion in July 2020. Within each plot, we installed a strip of 20 cm × 1 m and clipped all aboveground plant biomass. Samples were sorted to separate plant litter, monocots, forbs (including legumes), and mosses. All biomass was dried until constant weight at 60°C for at least 48hrs and weighed. We recorded plant community composition by visually estimating the percentage cover of each plant species using the decimal Londo scale (Londo 1976) in a 1 × 1 m permanent subplot within each plot. We used cover data to calculate plant species diversity as the Shannon index, and we used the biomass data to calculate the percentage of forb biomass relative to the total biomass of forbs and grasses.

We measured flower abundance by counting floral units in five 50 × 50 cm subplots per 5 × 5 m plot, for each flowering species, except wind-pollinated plants, such as grasses. Floral units are defined as approximately 1 cm² with at least one open flower (Carvalheiro et al. 2013). When the flowers are larger than 1 cm², the number of floral units corresponds to the flower size in cm², excluding petals. At the peak of the growing season, we also measured the amount of light reaching the soil surface by measuring light availability above and below the canopy on a sunny day between 11 am and 2 pm using a PAR meter (Li-COR). In each plot, we also recorded the area of bare soil by visually estimating the percentage cover by bare soil during the plant community assessments.

We quantified the abundance of different groups of insects in 2020, focussing on pollinators and herbivorous insects (specifically grasshoppers). For pollinators, we carried out flower visitation surveys between 11:30 hrs and 16:00 hrs on days with sunny to partly cloudy conditions, temperatures between 18 and 30°C, and a maximum wind speed of 5 Bft. Each 5 × 5 m plot was observed for 20 minutes by walking very slowly around the plots to minimise disturbance. A flower visitation event was defined as a moment when a pollinator made contact with the pollen or nectar source of the flower. For each visitation event, the plant species and the pollinator species were identified, either directly in the field or via photographs taken during the visitation event. A photo library was made for all observed pollinator species at the site, and identification was checked with taxonomic experts. Observed pollinators belonged to the following groups: Hymenoptera (bees and wasps), Diptera (hoverflies), and Lepidoptera (butterflies and moths). If an individual pollinator visited multiple plant species, only the visitation of the first plant species was registered. This method aligns with an unpublished Nutrient Network protocol where insects were caught upon visitation for identification in the laboratory. We used the pollinator observations to calculate Shannon diversity.

Grasshoppers (Orthoptera) were sampled on 12 August 2020 by sweep netting. Samples were collected between 10:00 hrs and 12:30 hrs. The temperature was circa 30°C and the windspeed was 1 Bft. In each plot, 50 sweeps were made with the net in 5 transects that were equally spaced and covering the whole 5 x 5 m plot. The transects were sampled from the outside of the plot towards the inside to prevent grasshoppers from escaping from the plot by the disturbance of the sweep netting. The sampling of one plot took approximately one minute. Samples were transferred to a plastic bag, which was placed in a -20°C freezer within 3 hours after sampling. For each sample, the number of adult individuals was counted and identified to species or genus level.