This should be a fairly straight forward but I still couldn't quite get it. Let $X$ be $n$ by $1$ random vector, $Y$ be another $n$ by $1$ and that $X$ and $Y$ are independent. What is $Var(X^tY)=:V_{XY}$?
I tried generalizing this formula from the univariate case $Var(xy)=E(x)^2Var(y)+E(y)^2Var(x)+Var(x)Var(y)$
except I hit a bit of a roadblock for the 3rd term above. That is, I get $$ V_{XY}=X^tV_YX+Y^tV_XY\\+E[(X^t-\mu_X^t)(Y-\mu_Y)(Y^t-\mu_Y^t)(X-\mu_X)] $$ So... what is it?
This answer supposes that $X^TY$ (where $X$ and $Y$ are $n\times 1$ vectors) is a $1\times 1$ vector or scalar $\sum_i X_iY_i$ and so we need to consider the variance of a single random variable that is this sum of products. Since $X$ and $Y$ are independent random vectors, we note that $X_1, Y_1$ are independent random variables as are $X_2, Y_2$. Also, $(X_1, X_2)$ is independent of $(Y_1, Y_2)$.
Since $\operatorname{var}(Z_1+Z_2) = \operatorname{var}(Z_1) + \operatorname{var}(Z_2) + 2\operatorname{cov}(Z_1,Z_2)$, we get that $$\begin{align} \operatorname{var}(X_1Y_1+X_2Y_2) &= \operatorname{var}(X_1Y_1)+ \operatorname{var}(X_2Y_2) + 2\operatorname{cov}(X_1Y_1, X_2Y_2)\\ \end{align}$$ The first two terms are the variances of the products of independent random variables, and the OP knows the formula for handling these. Turning to the covariance, we have that $$\begin{align} \operatorname{cov}(X_1Y_1, X_2Y_2) &= E[X_1Y_1X_2Y_2] - E[X_1Y_1]E[X_2Y_2]\\ &= E[X_1X_2]E[Y_1Y_2] - E[X_1]E[Y_1]E[X_2]E[Y_2]\\ &= \left(\operatorname{cov}(X_1, X_2)+E[X_1]E[X_2]\right)\left(\operatorname{cov}(Y_1, Y_2)+E[Y_1]E[Y_2]\right)\\ & \qquad \qquad - E[X_1]E[Y_1]E[X_2]E[Y_2]\\ \end{align}$$ which simplifies to $$ \operatorname{cov}(X_1Y_1, X_2Y_2) = \operatorname{cov}(X_1, X_2)E[Y_1]E[Y_2]+ \operatorname{cov}(Y_1, Y_2)E[X_1]E[X_2]\\ + \operatorname{cov}(X_1, X_2)\operatorname{cov}(Y_1, Y_2),\tag{1}$$ a result that is eerily similar to the result $$\operatorname{var}(X_iY_i) = \operatorname{var}(X_i)(E[Y_i])^2 + \operatorname{var}(Y_i)(E[X_i])^2 + \operatorname{var}(X_i)\operatorname{var}(Y_i) \tag{2}$$ quoted by the OP for independent random variables $X$ and $Y$
Thus we have $$\operatorname{var}\left(\sum_{i=1}^n X_iY_i\right) = \sum_{i=1}^n \operatorname{var}(X_iY_i) + 2\sum_{i=1}^{n-1}\sum_{j=i+1}^n \operatorname{cov}(X_iY_i, X_jY_j)$$ where each term in the first sum on the right is given by $(2)$ and each term in the double sum on the right is given by $(1)$. Plug and chug and enjoy!
Here's a matrix-notation version that might be more convenient to work with than the plug-and-chug required by Dilip's answer. $\DeclareMathOperator{\E}{\mathbb E} \DeclareMathOperator{\Var}{Var} \newcommand{\tp}{^\mathsf{T}} \DeclareMathOperator{\Tr}{Tr}$
Let $\mu_X := \E[X]$, $\Sigma_X := \Var[X]$, and similarly for $\mu_Y$, $\Sigma_Y$. Then, using trace rotation and independence of $X$ and $Y$, we get \begin{align} \Var[X\tp Y] &= \E[(X\tp Y)^2] - \E[X\tp Y]^2 \\&= \E\left[ X\tp Y Y\tp X\right] - \left(\E[X]\tp \E[Y]\right)^2 \\&= \E\left[\Tr\left( X X\tp Y Y\tp \right)\right] - \left(\mu_X\tp \mu_Y\right)^2 \\&= \Tr\left( \E[X X\tp] \E[Y Y\tp] \right) - \left(\mu_X\tp \mu_Y\right)^2 \\&= \Tr\left( \left( \mu_X \mu_X\tp + \Sigma_X \right) \left( \mu_Y \mu_Y\tp + \Sigma_Y \right) \right) - \left(\mu_X\tp \mu_Y\right)^2 \\&= \Tr\left( \mu_X \mu_X\tp \mu_Y \mu_Y\tp \right) + \Tr\left( \mu_X \mu_X\tp \Sigma_Y \right) + \Tr\left( \Sigma_X \mu_Y \mu_Y\tp \right) + \Tr\left( \Sigma_X \Sigma_Y \right) - \left(\mu_X\tp \mu_Y\right)^2 \\&= \left(\mu_X\tp \mu_Y\right)^2 + \Tr\left( \mu_X\tp \Sigma_Y \mu_X \right) + \Tr\left( \mu_Y\tp \Sigma_X \mu_Y \right) + \Tr\left( \Sigma_X \Sigma_Y \right) - \left(\mu_X\tp \mu_Y\right)^2 \\&= \mu_X\tp \Sigma_Y \mu_X + \mu_Y\tp \Sigma_X \mu_Y + \Tr(\Sigma_X \Sigma_Y) .\end{align}
