[proofplan]
Represent the Wishart matrix as a sum of independent rank-one matrices $X_rX_r^\top$. The expectation follows by taking expectations entrywise and using the covariance matrix of a centred Gaussian vector. For the covariance, expand each entry $w_{ij}$ as a sum over observations; independence removes all cross-observation covariance terms, and the remaining single-observation fourth moment is evaluated by Isserlis' formula. The variance statement is the special case $(k,l)=(i,j)$ of the covariance formula.
[/proofplan]
[step:Expand the Wishart entries as sums of Gaussian products]
For each $r \in \{1,\dots,n\}$, the outer product $X_rX_r^\top$ is the $p \times p$ random matrix whose $(i,j)$ entry is $X_{ri}X_{rj}$. Hence, for every $i,j \in \{1,\dots,p\}$,
\begin{align*}
w_{ij}
=
\sum_{r=1}^n X_{ri}X_{rj}.
\end{align*}
Since $X_r \sim \mathcal{N}_p(0,\Sigma)$, its coordinate means and covariances satisfy
\begin{align*}
\mathbb{E}[X_{ri}] &= 0,\\
\mathbb{E}[X_{ri}X_{rj}] &= \sigma_{ij}.
\end{align*}
[/step]
[step:Compute the expectation entry by entry]
Fix $i,j \in \{1,\dots,p\}$. By linearity of expectation applied to the finite sum defining $w_{ij}$,
\begin{align*}
\mathbb{E}[w_{ij}]
&=
\mathbb{E}\left[\sum_{r=1}^n X_{ri}X_{rj}\right] \\
&=
\sum_{r=1}^n \mathbb{E}[X_{ri}X_{rj}] \\
&=
\sum_{r=1}^n \sigma_{ij} \\
&=
n\sigma_{ij}.
\end{align*}
Since this holds for every matrix entry, the matrix identity is
\begin{align*}
\mathbb{E}[W] = n\Sigma.
\end{align*}
[/step]
[step:Reduce the covariance to one Gaussian observation]
Fix $i,j,k,l \in \{1,\dots,p\}$. Define real-valued random variables
\begin{align*}
A_r &: \Omega \to \mathbb{R}, & A_r(\omega) &= X_{ri}(\omega)X_{rj}(\omega),\\
B_r &: \Omega \to \mathbb{R}, & B_r(\omega) &= X_{rk}(\omega)X_{rl}(\omega),
\end{align*}
for $r \in \{1,\dots,n\}$. Then
\begin{align*}
w_{ij} = \sum_{r=1}^n A_r,
\qquad
w_{kl} = \sum_{r=1}^n B_r.
\end{align*}
By bilinearity of covariance over finite sums,
\begin{align*}
\operatorname{Cov}(w_{ij},w_{kl})
=
\sum_{r=1}^n \sum_{s=1}^n \operatorname{Cov}(A_r,B_s).
\end{align*}
If $r \neq s$, then $X_r$ and $X_s$ are independent, so $A_r$ and $B_s$ are independent. Therefore
\begin{align*}
\operatorname{Cov}(A_r,B_s)=0
\qquad
\text{for } r \neq s.
\end{align*}
Thus only the diagonal terms remain:
\begin{align*}
\operatorname{Cov}(w_{ij},w_{kl})
=
\sum_{r=1}^n \operatorname{Cov}(A_r,B_r).
\end{align*}
Because the random vectors $X_1,\dots,X_n$ have the same distribution, the summands are equal, and hence
\begin{align*}
\operatorname{Cov}(w_{ij},w_{kl})
=
n\,\operatorname{Cov}(X_{1i}X_{1j},X_{1k}X_{1l}).
\end{align*}
[guided]
The entry $w_{ij}$ is a sum over the $n$ independent observations, and the same is true of $w_{kl}$. To make the covariance expansion explicit, define
\begin{align*}
A_r &: \Omega \to \mathbb{R}, & A_r(\omega) &= X_{ri}(\omega)X_{rj}(\omega),\\
B_r &: \Omega \to \mathbb{R}, & B_r(\omega) &= X_{rk}(\omega)X_{rl}(\omega).
\end{align*}
Then
\begin{align*}
w_{ij} = \sum_{r=1}^n A_r,
\qquad
w_{kl} = \sum_{r=1}^n B_r.
\end{align*}
Using bilinearity of covariance for finite sums gives
\begin{align*}
\operatorname{Cov}(w_{ij},w_{kl})
=
\sum_{r=1}^n \sum_{s=1}^n \operatorname{Cov}(A_r,B_s).
\end{align*}
The purpose of introducing $A_r$ and $B_s$ is to separate terms coming from the same observation from terms coming from different observations. If $r \neq s$, then $A_r$ is a measurable function of $X_r$ and $B_s$ is a measurable function of $X_s$. Since $X_r$ and $X_s$ are independent, $A_r$ and $B_s$ are independent, and hence
\begin{align*}
\operatorname{Cov}(A_r,B_s)
=
\mathbb{E}[A_rB_s]-\mathbb{E}[A_r]\mathbb{E}[B_s]
=
0.
\end{align*}
Therefore all off-diagonal terms in the double sum vanish:
\begin{align*}
\operatorname{Cov}(w_{ij},w_{kl})
=
\sum_{r=1}^n \operatorname{Cov}(A_r,B_r).
\end{align*}
Since $X_1,\dots,X_n$ are identically distributed, each covariance in the last sum has the same value. Thus
\begin{align*}
\operatorname{Cov}(w_{ij},w_{kl})
=
n\,\operatorname{Cov}(X_{1i}X_{1j},X_{1k}X_{1l}).
\end{align*}
[/guided]
[/step]
[step:Evaluate the single-observation fourth moment by Isserlis' formula]
Let $Y : \Omega \to \mathbb{R}^p$ be the random vector $Y := X_1$, and write $Y = (Y_1,\dots,Y_p)$. Then $Y \sim \mathcal{N}_p(0,\Sigma)$, so
\begin{align*}
\mathbb{E}[Y_a] = 0,
\qquad
\mathbb{E}[Y_aY_b] = \sigma_{ab}
\end{align*}
for all $a,b \in \{1,\dots,p\}$.
By Isserlis' formula for centred multivariate Gaussian fourth moments (citing a result not yet in the wiki: Isserlis' formula), applied to the four coordinates $Y_i,Y_j,Y_k,Y_l$,
\begin{align*}
\mathbb{E}[Y_iY_jY_kY_l]
=
\mathbb{E}[Y_iY_j]\mathbb{E}[Y_kY_l]
+
\mathbb{E}[Y_iY_k]\mathbb{E}[Y_jY_l]
+
\mathbb{E}[Y_iY_l]\mathbb{E}[Y_jY_k].
\end{align*}
Substituting the covariance entries of $Y$ gives
\begin{align*}
\mathbb{E}[Y_iY_jY_kY_l]
=
\sigma_{ij}\sigma_{kl}
+
\sigma_{ik}\sigma_{jl}
+
\sigma_{il}\sigma_{jk}.
\end{align*}
Therefore
\begin{align*}
\operatorname{Cov}(Y_iY_j,Y_kY_l)
&=
\mathbb{E}[Y_iY_jY_kY_l]
-
\mathbb{E}[Y_iY_j]\mathbb{E}[Y_kY_l] \\
&=
\left(\sigma_{ij}\sigma_{kl}
+
\sigma_{ik}\sigma_{jl}
+
\sigma_{il}\sigma_{jk}\right)
-
\sigma_{ij}\sigma_{kl} \\
&=
\sigma_{ik}\sigma_{jl}
+
\sigma_{il}\sigma_{jk}.
\end{align*}
[guided]
We now need the covariance of two quadratic coordinate products from one centred Gaussian vector. Set $Y := X_1$, so $Y : \Omega \to \mathbb{R}^p$ is a random vector with distribution $\mathcal{N}_p(0,\Sigma)$, and write $Y=(Y_1,\dots,Y_p)$. The covariance matrix condition means
\begin{align*}
\mathbb{E}[Y_a] = 0,
\qquad
\mathbb{E}[Y_aY_b] = \sigma_{ab}
\end{align*}
for every $a,b \in \{1,\dots,p\}$.
The relevant fourth moment is computed by Isserlis' formula for centred Gaussian random variables (citing a result not yet in the wiki: Isserlis' formula). Applied to the four centred Gaussian coordinates $Y_i,Y_j,Y_k,Y_l$, it gives
\begin{align*}
\mathbb{E}[Y_iY_jY_kY_l]
=
\mathbb{E}[Y_iY_j]\mathbb{E}[Y_kY_l]
+
\mathbb{E}[Y_iY_k]\mathbb{E}[Y_jY_l]
+
\mathbb{E}[Y_iY_l]\mathbb{E}[Y_jY_k].
\end{align*}
Each expectation on the right is an entry of the covariance matrix $\Sigma$, so
\begin{align*}
\mathbb{E}[Y_iY_jY_kY_l]
=
\sigma_{ij}\sigma_{kl}
+
\sigma_{ik}\sigma_{jl}
+
\sigma_{il}\sigma_{jk}.
\end{align*}
Now subtract the product of the means of the two quadratic factors:
\begin{align*}
\operatorname{Cov}(Y_iY_j,Y_kY_l)
&=
\mathbb{E}[Y_iY_jY_kY_l]
-
\mathbb{E}[Y_iY_j]\mathbb{E}[Y_kY_l] \\
&=
\left(\sigma_{ij}\sigma_{kl}
+
\sigma_{ik}\sigma_{jl}
+
\sigma_{il}\sigma_{jk}\right)
-
\sigma_{ij}\sigma_{kl} \\
&=
\sigma_{ik}\sigma_{jl}
+
\sigma_{il}\sigma_{jk}.
\end{align*}
The cancellation of $\sigma_{ij}\sigma_{kl}$ is exactly the subtraction built into covariance.
[/guided]
[/step]
[step:Conclude the covariance and variance formulas]
Combining the reduction to one observation with the single-observation covariance calculation gives
\begin{align*}
\operatorname{Cov}(w_{ij},w_{kl})
&=
n\,\operatorname{Cov}(X_{1i}X_{1j},X_{1k}X_{1l}) \\
&=
n(\sigma_{ik}\sigma_{jl}+\sigma_{il}\sigma_{jk}).
\end{align*}
This proves the asserted entrywise covariance formula.
Taking $k=i$ and $l=j$ gives
\begin{align*}
\operatorname{Var}(w_{ij})
&=
\operatorname{Cov}(w_{ij},w_{ij}) \\
&=
n(\sigma_{ii}\sigma_{jj}+\sigma_{ij}\sigma_{ji}).
\end{align*}
Since $\Sigma$ is symmetric, $\sigma_{ji}=\sigma_{ij}$, and hence
\begin{align*}
\operatorname{Var}(w_{ij})
=
n(\sigma_{ii}\sigma_{jj}+\sigma_{ij}^2).
\end{align*}
Together with $\mathbb{E}[W]=n\Sigma$, this completes the proof.
[/step]
