This should readily follow if you look at the expression for the moment generating function (M.G.F.) for a multivariate normal density. M.G.F. looks like:
$E_X(e^{t'X}) = e^{\mu't + 0.5t'\Sigma t}$ for any real $n$-vector $t$, where $X = (X_1, X_2,...,X_n)$ is multivariate normal $N_n(\mu,\Sigma)$ random variate [simple to derive this!]
You can see from the expression that the exponent $\mu't + 0.5t'\Sigma t$ is a scalar and moreover, if you want to find the M.G.F. of a $p$-subset ($p \ne 0$) of the random vector $X$, say, {${X_{i_1}, X_{i_2}, ..., X_{i_p}}$}, where, {$i_1,i_2,..,i_p$} is a permutation of {$1,2,..,n$}, then you just need to plug in $t_j=0$ for $j$ $\epsilon$ {$1,2,..,n$}-{$i_1,i_2,..,i_p$} in the expression of M.G.F. The simplified form reduces to a known M.G.F. of a uni/multi-variate normal and that M.G.F. uniquely determines a distribution confirms the nice property of the marginal density of multivariate normal distribution.
You need to perform a linear change of variables to do this. Let $H$ be a solution to $H^TH=\Sigma^{-1}$ (it exists because $\Sigma$ is positive definite). We will consider the random vector $y=H(x-\mu)$.
We have $y^Ty=(x-\mu)^TH^TH(x-\mu)=(x-\mu)^T\Sigma^{-1}(x-\mu)$ and $|H||\Sigma|^{1/2}=1$. So the pdf of $y$ is given by
$p(y)=\frac1{|H|}p(x)=\frac{1}{|H|}\frac{1}{(2\pi)^{n/2}|\Sigma|^{1/2}}\exp\left(-\frac12(x-\mu)^T\Sigma^{-1}(x-\mu)\right)= \frac{1}{(2\pi)^{n/2}}\exp\left(-\frac12y^Ty\right)$
This means that the entries of $y$ are independent standard normal variables, so $\mathrm{Cov}[y]=I$.
In general we have $\mathrm{Cov}[Ax+b]=A\mathrm{Cov}[X]A^T$ (this is the multivariate version of $\mathrm{Var}(ax+b)=a^2\mathrm{Var}(x)$, and can be shown straightforwardly by considering the definition $\mathrm{Cov}[X]=\mathbb{E}[XX^T]-\mathbb{E}[X]\mathbb{E}[X]^T$).
So in our case $I=\mathrm{Cov}[y]=\mathrm{Cov}[H(x-\mu)]=H\mathrm{Cov}[x]H^T$, and $\mathrm{Cov}[x]=H^{-1}(H^T)^{-1}=(H^TH)^{-1}=\Sigma$.
Best Answer
It is easiest to look at this from the characteristic function point of view.
$\phi(t_1, t_2, \ldots, t_n) = \mathbb{E}( \exp(i \langle \vec{t}, \vec{x} \rangle)) = \exp( i \langle \vec{t}, \vec{\mu}\rangle - \frac{1}{2} \langle \vec{t}, \Sigma \vec{t}\rangle)$.
Characteristic function of the marginal distribution is the characteristic function of the original multi-normal with the marginalized $t$-s being zero.
Clearly setting $t_i = 0$ is equivalent to dropping out $i$-th variable.