I think this is not quite what you had in mind, but for the sake of completeness, let's start with some basics. Most quadrature formulas such as Newton-Cotes and Gauss are based on the idea that in order to evaluate the integral of a function approximately, you can approximate the function by, e.g., a polynomial that you can then integrate exactly:
$$ \int_a^b f(x) \,dx \approx \int_a^b \sum_j c_j p_j(x)\,dx = \sum_j c_j\int_a^b p_j(x) \,dx. $$
Newton-Cotes and Gauss are based on Lagrange interpolation, meaning you interpolate the given function using its values on a set of nodes $x_j$ (which are spaced uniformly for Newton-Cotes and chosen optimally in a certain sense for Gauss). In this case, $c_j = f(x_j)$, and the integrals over the polynomial nodal basis functions $p_j$ are exactly the quadrature weights.
The same approach works with Hermite interpolation, i.e., interpolation using the values of a function and its derivatives up to a certain order on a set of nodes. In the case of function and first derivative values only, you have
$$ \int_a^b f(x) \,dx \approx \int_a^b \sum_j f(x_j) p_j(x) + f'(x_j) q_j(x)\,dx = \sum_j f(x_j) w_j+f'(x_j) \bar w_j. $$
(There is a Matlab implementation of this, if you'd like to see how it works.)
This is related to a variant of Gauss quadrature called Gauss-Legendre quadrature, where the nodes are chosen precisely to make the weights $\bar w_j$ vanish (which is another explanation for the fact that Gauss quadrature with $N$ nodes is exact of order $2N-1$). I think this at least partially answers your question in the second paragraph.
For this reason, Gauss quadrature is usually used instead of Hermite interpolation, since you get the same order with the same number of points, but do not need derivative information.
For multidimensional quadrature, you face the problem that the number of derivatives (including mixed derivatives) you need to evaluate grows very quickly as the order increases.
Coming back to your question: A straightforward way of exploiting derivative information would be to use a subdivision of your integration domain and use a separate quadrature for every division. If you know that the derivatives of your function are large in some part of the domain, you'd use either smaller domains (in effect, a summed quadrature formula) or higher quadrature order. This is related to h- and p-adaptivity, respectively, in finite element methods.