Quasi-interpolation is an important tool, used both in theory and in practice, for the approximation of smooth functions from univariate or multivariate spaces which contain Π m = Π m ( R d ) , the d -variate polynomials of degree ≤ m . In particular, the reproduction of Π m leads to an approximation order of m + 1 . Prominent examples include Lagrange and Bernstein type approximations by polynomials, the orthogonal projection onto Π m for some inner product, finite element methods of precision m , and multivariate spline approximations based on macroelements or the translates of a single spline. For such a quasi-interpolation operator L which reproduces Π m ( R d ) and any r ≥ 0 , we give an explicit construction of a quasi-interpolant R m r + m L = L + A which reproduces Π m + r , together with an integral error formula which involves only the ( m + r + 1 ) th derivative of the function approximated. The operator R m m + r L is defined on functions with r additional orders of smoothness than those on which L is defined. This very general construction holds in all dimensions d . A number of representative examples are considered.