Abstract It is well known that Gaussian cubature rules are related to multivariate orthogonal polynomials. The cubature rules found in the literature use common zeroes of some linearly independent set of products of basically univariate polynomials. We show how a new family of multivariate orthogonal polynomials, so-called spherical orthogonal polynomials, leads to symbolic–numeric Gaussian cubature rules in a very natural way. They can be used for the integration of multivariate functions that in addition may depend on a vector of parameters and they are exact for multivariate parameterized polynomials. Purely numeric Gaussian cubature rules for the exact integration of multivariate polynomials can also be obtained. We illustrate their use for the symbolic–numeric solution of the partial differential equations satisfied by the Appell function F 2 , which arises frequently in various physical and chemical applications. The advantage of a symbolic–numeric formula over a purely numeric one is that one obtains a continuous extension, in terms of the parameters, of the numeric solution. The number of symbolic–numeric nodes in our Gaussian cubature rules is minimal, namely m for the exact integration of a polynomial of homogeneous degree 2 m − 1 . In Section 1 we describe how the symbolic–numeric rules are constructed, in any dimension and for any order. In Sections 2, 3 and 4 we explicit them on different domains and for different weight functions. An illustration of the new formulas is given in Section 5 and we show in Section 6 how numeric cubature rules can be derived for the exact integration of multivariate polynomials. From Section 7 it is clear that there is a connection between our symbolic–numeric cubature rules and numeric cubature formulae with a minimal (or small) number of nodes.