Although the original question is a little too narrow, I'm going to interpret it as asking about general multidimensional integrals in which Gaussians appear together with arbitrary factors in the integrand:

$$\iiint f(\,\vec{r}\,)\,\exp(\,-\frac{1}{2}\vec{r}^T\Sigma^{-1}\vec{r}\,)\,dx\,dy\,dz$$

Here $\Sigma^{-1}$ could be a symmetric positive definite matrix, or in the simplest case just $\sigma^2$ times the identity matrix. The dimension can be anything, but here I'll just specialize to three.

Then you may be interested in this observation which speeds up such integrals tremendously in comparison with straightforward multiple integration:

**Edit: An integrand that doesn't factor**

Start by defining the Gaussian in three dimensions:

    f = PDF[MultinormalDistribution[{0, 0, 
        0}, \[Sigma] IdentityMatrix[3]], {x, y, z}]

> $\frac{e^{\frac{1}{2}
   \left(-\frac{x^2}{\sigma
   }-\frac{y^2}{\sigma
   }-\frac{z^2}{\sigma }\right)}}{2
   \sqrt{2} \pi ^{3/2} \sqrt{\sigma
   ^3}}$

Try the following as a triple integral:

    Timing[Integrate[
      Integrate[Integrate[
        (x^2 + y^3)/(1 + z^2)
          f, {x, -Infinity, Infinity}], {y, -Infinity, \
    Infinity}], {z, -Infinity, Infinity}]
     ]

> $\left\{1.48519,\text{ConditionalExpression}\left[\frac{\sqrt{\frac{
   \pi }{2}} e^{\frac{1}{2 \sigma
   }} \sqrt{\sigma ^3}
   \text{erfc}\left(\frac{1}{\sqrt{
   2} \sqrt{\sigma
   }}\right)}{\sigma },\Re(\sigma
   )>0\right]\right\}$

Now do the same using `Expectation`, keeping in mind the different normalization of the result:

    Timing[Expectation[(x^2 + y^3)/(
      1 + z^2), {x, y, z} \[Distributed] 
       MultinormalDistribution[{0, 0, 0}, \[Sigma] IdentityMatrix[3]]]
     ]

> $\left\{0.775868,\sqrt{\frac{\pi
   }{2}} e^{\frac{1}{2 \sigma }}
   \sqrt{\sigma }
   \text{erfc}\left(\frac{1}{\sqrt{
   2} \sqrt{\sigma
   }}\right)\right\}$

So the result is obtained about twice as fast using `Expectation`. 

**Fourier transforms as a special case**

My initial example was the Fourier transform of a Gaussian, 

$$\iiint \exp(- i\, \vec{k}\cdot \vec{r}\,)\,\exp(\,-\frac{1}{2}\vec{r}^T\Sigma^{-1}\vec{r}\,)\,dx\,dy\,dz$$

where the speed advantage of `Expectation` compared to multiple `Integrate` with infinite integration boundaries is more significant than above:

    Expectation[
     Exp[-I {kx, ky, kz}.{x, y, z}], {x, y, z} \[Distributed] 
      MultinormalDistribution[{0, 0, 0}, σ IdentityMatrix[3]]]

> $e^{-\frac{1}{2} \sigma 
   \left(\text{kx}^2+\text{ky}^2+\text{kz}^2\right)}$

This works for *symbolic* integrands (note I didn't specify a value for $\sigma$ or the components of $\vec{k}$), and you don't have to specify the list of integration boundaries explicitly either. The timing for the above is about 7 times faster than doing the calculation using three nested `Integrate`s, where I also have to add assumptions for the parameters (specifically $\sigma>0$, $\vec{k}$ real).

Although my initial claim is correct that this is faster than the straightforward triple integral, one can do a lot better in this special case by just using `FourierTransform` directly:

    Timing[FourierTransform[f, {x, y, z}, {kx, ky, kz}]]

> $\left\{0.007332,\frac{e^{-\frac{1}{2} \sigma 
   \left(\text{kx}^2+\text{ky}^2+\text{kz}^2\right)}
   }{2 \sqrt{2} \pi ^{3/2} \left(\frac{1}{\sigma
   }\right)^{3/2} \sqrt{\sigma ^3}}\right\}$

Interestingly, though, we can *even* beat this really fast result using the multinormal distribution again:

    Timing[CharacteristicFunction[
      MultinormalDistribution[{0, 0, 0}, \[Sigma] IdentityMatrix[3]], {kx,
        ky, kz}]
     ]

> $\left\{0.000185,e^{\frac{1}{2} \left(\text{kx}^2
   (-\sigma )-\text{ky}^2 \sigma -\text{kz}^2 \sigma
   \right)}\right\}$

**Conclusion**

For Gaussian integrals over all space (or momentum space, as in the question), the  approach using `MultinormalDistribution` is complementary to whuber's solution: general Gaussian integrals can be evaluated by using `Expectation` and similar tools for probability distributions, such as `CharacteristicFunction`.