As the passage from step (c) to step (d) is non-trivial, let me add a new answer to this interesting question here.
(b) Let $A = O'\Lambda O$ be the spectral decomposition of $A$, where $O$ is an order $n$ orthogonal matrix and $\Lambda = \operatorname{diag}(\lambda_1, \ldots, \lambda_n)$. Let $z = Oy$, then $z \sim N(0, I_n)$ as $y \sim N(0, I_n)$ and $O$ is orthogonal. It then follows $y'Ay = y'O'\Lambda Oy = z'\Lambda z = \sum\limits_{i = 1}^n\lambda_i z_i^2$ and the independence of $z_1^2, \ldots, z_n^2$ that
\begin{align}
\operatorname{Var}(y'Ay) = \sum_{i = 1}^n\lambda_i^2\operatorname{Var}(z_i^2) =
\sum_{i = 1}^n\lambda_i^2(E[z_i^4] - (E[z_i^2])^2) = 2\sum_{i = 1}^n\lambda_i^2 =
2\operatorname{tr}(A^2). \tag{1}
\end{align}
In $(1)$, we used that $E[z_i^4] = 3$ and $E[z_i^2] = 1$ if $z_i \sim N(0, 1)$.
(c) If $y \sim N(0, \Sigma)$, we can write $y = \Sigma^{1/2}z$, where $z \sim N(0, I_n)$, it then follows by $(1)$ and the symmetry of $(\Sigma^{1/2})'A\Sigma^{1/2}$ that
\begin{align}
\operatorname{Var}(y'Ay) &= \operatorname{Var}(z'(\Sigma^{1/2})'A\Sigma^{1/2}z) =
2\operatorname{tr}((\Sigma^{1/2}A\Sigma^{1/2})^2) \\
&= 2\operatorname{tr}(\Sigma^{1/2}A\Sigma A\Sigma^{1/2})
= 2\operatorname{tr}(A\Sigma A\Sigma) = 2\operatorname{tr}((A\Sigma)^2). \tag{2}
\end{align}
In $(2)$, we used that $(\Sigma^{1/2})' = \Sigma^{1/2}$ and $\operatorname{tr}(M_1M_2) = \operatorname{tr}(M_2M_1)$ for order $n$ matrices $M_1, M_2$.
(d) If $y \sim N(\mu, \Sigma)$, we can write $y = \mu + z$, where $z \sim N(0, \Sigma)$, whence
\begin{align}
y'Ay = (\mu + z)'A(\mu + z) = \mu'A\mu + 2\mu'Az + z'Az.
\end{align}
Therefore,
\begin{align}
\operatorname{Var}(y'Ay) &= \operatorname{Var}(z'Az + 2\mu'Az + \mu'A\mu) = \operatorname{Var}(z'Az + 2\mu'Az) \\
&=\operatorname{Var}(z'Az) + 4\operatorname{Var}(\mu'Az) + 2\operatorname{Cov}(z'Az,
2\mu'Az). \tag{3}
\end{align}
By $(2)$,
\begin{align}
\operatorname{Var}(z'Az) = 2\operatorname{tr}((A\Sigma)^2). \tag{4}
\end{align}
In addition,
\begin{align}
\operatorname{Var}(\mu'Az) = \mu'A\operatorname{Var}(z)A\mu = \mu'A\Sigma A\mu.
\tag{5}
\end{align}
To evaluate $\operatorname{Cov}(z'Az, 2\mu'Az) = 2E[z'Az\mu'Az]$, note that the expansion of $z'Az\mu'Az$ consists of only terms of forms (regardless constant coefficients) $z_i^3, 1 \leq i \leq n$ and $z_i^2z_j, 1 \leq i \neq j \leq n$, while for $z \sim N(0, \Sigma)$, it is easy to verify that
\begin{align}
& E[z_i^3] = 0, \quad 1 \leq i \leq n, \\
& E[z_i^2z_j] = 0, \quad 1 \leq i \neq j \leq n.
\end{align}
Therefore,
\begin{align}
\operatorname{Cov}(z'Az, 2\mu'Az) = 0. \tag{6}
\end{align}
Substituting $(4), (5), (6)$ into $(3)$, we obtain
\begin{align}
\operatorname{Var}(y'Ay) = 2\operatorname{tr}((A\Sigma)^2) + 4\mu'A\Sigma A\mu.
\end{align}
This completes the proof.