Suppose $\vec x \sim N(\vec \mu, \Sigma)$ is multivariate normal. I want to see that $\mathrm{Cov}(\vec x^TA\vec x,\vec x^TB\vec x) = 2 \mathrm{Tr}(A \Sigma B \Sigma) + 4 \vec \mu^TA \Sigma B \vec \mu$. I have been searching the internet for a while, and found multiple sources confirming it, but no proof, for example here: $Var(Q)=2\ tr(A\Sigma A\Sigma)+4\mu^TA\Sigma A\mu$

Do you know how to show it?