Abstract In this paper, we study the Sobolev regularity of solutions to nonlinear second order elliptic equations with super-linear first-order terms on Riemannian manifolds, complemented with Neumann boundary conditions, when the source term of the equation belongs to a Lebesgue space, under various integrability regimes. Our method is based on an integral refinement of the Bochner identity, and leads to “semilinear Calderón–Zygmund” type results. Applications to the problem of smoothness of solutions to Mean Field Games systems with Neumann boundary conditions posed on convex domains of the Euclidean space will also be discussed.