If we have a path integral involving many fields,
$$Z = \int \mathcal D \phi_1 \cdots \mathcal D \phi_n \exp(-S[\phi_1,\ldots, \phi_n]),$$
and $\phi_n$ occurs only quadratically-- i.e. the $\mathcal D \phi_n$ integral is Gaussian-- I've been told that integrating over $\phi_n$ is equivalent to solving for $\phi_n$'s equation of motion
$$\phi_n= f(\phi_1,\ldots, \phi_{n-1})$$
using Euler-Lagrange and plugging in. Up to normalization. Can one show in general why this is true?
The gaussian integral $$ \int dx\,e^{-\frac12 a x^2 + bx + c} = \sqrt{\frac{2\pi}{a}}\, e^{c+b^2/(2a)}\,, $$ is similar to its path integral counterpart, which is $$ \int \mathcal{D}\phi\,e^{-\frac12\phi \cdot A\cdot \phi + \phi\cdot b + C} \propto \exp\left(C + \frac{1}{2} \, b\cdot A^{-1}\cdot b\right)\,. $$ By the dot I mean $a\cdot b \equiv \int a(x)\, b(x)$, $a\cdot B \cdot c \equiv \int a(x)\, B(x,y)\, c(y)$. Moreover $A^{-1}$ satisfies $$ \int A(x,y)\cdot A^{-1}(y,z) = \delta(x-z)\,. $$
The equations of motion for $\phi$ are $$ -A\cdot \phi + b = 0\qquad \Longrightarrow\qquad\phi = A^{-1}\cdot b\,. $$ Replacing this on the action yields the same result $$ -\frac12\phi \cdot A\cdot \phi + \phi\cdot b + C \quad\to\quad -\frac12\,b \cdot A^{-1} \cdot A\cdot A^{-1}\cdot b + b\cdot A^{-1}\cdot b + C = \frac12\,b\cdot A^{-1}\cdot b + C\,. $$ If the dot notation is confusing I suggest to expand everything in integrals. The operator $A$ usually is just $(\square_x + m^2) \delta(x-y)$ and $A^{-1}$ is $G_F(x-y)$, the Feynman propagator.
Edit: As a comment pointed out, this does not take into account the $(\det A)^{-1/2}$. If $A$ is a constant operator, this does not pose any problem in perturbative computations because we only need the partition function modulo overall factors.
On the other hand, if $A$ is a function of the remaining fields $A(\phi_1,\ldots,\phi_{n-1})$, it will not pass through the subsequent integrals. The way this is normally handled is by exponentiating it as $$ (\det A)^{-1/2} = e^{- \frac12 \mathrm{Tr}\log A}\,, $$ (with a suitable regularization procedure) and this typically yields a non-local action
Gaussian integration is a particularly simple case of the WKB expansion, cf. e.g. this Phys.SE post. Of course, the caveat is that the saddle point may be complex-valued. In other words, in the 1D case, the saddle point may lie in the complex plane, and one has to show that one can close the integration contour between the real axis and the line of steepest descent through the saddle point. Some of these issues are addressed in e.g. this & this related posts.
