I don't understand the subscript $n$ notation, however, in the least squares problem that is given by:
\begin{equation}
{\bf{y}}={\bf{H}}{\theta}+\bf{n},
\end{equation}
where ${\bf{n}}\sim\mathcal{N}(\bf{0}, \sigma^2I_N)$ is a zero mean additive white Gaussian noise and $I_N$ is the $N \times N$ identity matrix, the maximum likelihood and the least squares estimators are equivalent and are efficient estimators (i.e. unbiased and attains the Cramer Rao bound) and are given by
\begin{equation}
\widehat \theta_{ML}=\left ( \bf{H^TH} \right )^{-1} \bf{H^Ty} .
\end{equation}
Now, let's calculate the first two moments of the estimator:
\begin{equation}
E\left \{ \widehat \theta_{ML} \right \} = E\left \{ \left ( \bf{H^TH} \right )^{-1} \bf{H^Ty} \right \} = \\
\left ( \bf{H^TH} \right )^{-1} \bf{H^T} \underset{=\bf{H}\theta}{\underbrace{E\left \{ \bf{y}\right \}}} = \theta,
\end{equation}
where $\theta$ is the true parameter value.
Thus, the estimator is unbiased and its Mean Square Error(MSE) will equal its covariance.
Now for the covariance:
Remember that if $\bf{z}=\bf{Ax}$ then $cov\left ( \bf{z} \right ) = \bf{A} cov\left ( x \right ) \bf{A^T}$.
Since the mean doesn't change the covariance in our problem $\bf{y} \sim \mathcal{N}\left ( \bf{H\theta}, \sigma^2I_N \right )$, thus,
\begin{equation}
cov\left ( \widehat \theta_{ML} \right )=\left ( \bf{H^TH} \right )^{-1} \bf{H^T} cov\left ( \bf{n} \right ) \bf{H}\left ( \bf{H^TH} \right )^{-1} = \\
= \left ( \bf{H^TH} \right )^{-1} \bf{H^T} \sigma^2 \bf{I_N} \bf{H}\left ( \bf{H^TH} \right )^{-1} = \\
= \sigma^2 \left ( \bf{H^TH} \right )^{-1} \bf{H^T} \bf{H}\left ( \bf{H^TH} \right )^{-1} = \sigma^2 \left ( \bf{H^TH} \right )^{-1}.
\end{equation}
Note that $\bf{H^TH}$ is symmetric so it equals its transpose. Same is correct for its inverse.
I believe that this answers your question completely.