The relationship between biodiversity and ecosystem functions (BEFs) has attracted great interest. Studies on BEF have so far focused on the average trend of ecosystem function as species diversity increases. A tantalizing but rarely addressed question is why large variations in ecosystem functions are often observed across systems with similar species diversity, likely obscuring observed BEFs. Here we use a multi-trophic food web model in combination with empirical data to examine the relationships between species richness and the variation in ecosystem functions (VEFs) including biomass, metabolism, decomposition, and primary and secondary production. We then probe the mechanisms underlying these relationships, focusing on the role of trophic interactions. While our results reinforce the previously documented positive BEF relationships, we found that ecosystem functions exhibit significant variation within each level of species richness and the magnitude of this variation displays a hump-shaped relationship with species richness. Our analyses demonstrate that VEFs is reduced when consumer diversity increases through elevated nonlinearity in trophic interactions, and/or when the diversity of basal species such as producers and decomposers decreases. This explanation is supported by a 34-year empirical food web time series from the Gulf of Riga ecosystem. Our work suggests that biodiversity loss may not only result in ecosystem function decline, but also reduce the predictability of functions by generating greater function variability among ecosystems. It thus helps to reconcile the debate on the generality of positive BEF relationships and to disentangle the drivers of ecosystem stability. The role of trophic interactions and the variation in their strengths mediated by functional responses in shaping ecosystem function variation warrants further investigations and better incorporation into biodiversity-ecosystem functioning research.