Consider the following multivariate (response) linear model:
\[\bm y_i = \bm x_i^T \bm\beta + \bm \epsilon_i,\]where $\bm y_i$ is $(1 \times r)$, $\bm x_i$ is $(1 \times p)$, $\bm \beta$ is $(p \times r)$, $\bm \epsilon_i$ is $(1\times r)$, $\text{Var}(\bm\epsilon_i) = \bm \Sigma$, and $\text{Cov}(\bm\epsilon_i, \bm\epsilon_j) = \bm 0$.
For $N$ (independent) data points, we can write this model as:
\[\begin{aligned} \vec\p{\bm Y} &= \vec\p{\bm X \bm\beta} + \vec\p{\bm\Xi} \end{aligned}\]where $\vec\p{\cdot}$ denotes column-stacking a matrix, $Y$ ($N\times r$) is constructed by stacking $N$ rows of $\bm y_i$, $\bm X~(N\times p)$ is constructed by stacking $N$ rows of $\bm x_i$, and $\bm\Xi$ is constructed by stacking $N$ rows of $\bm\epsilon_i$. Thus,
\[\begin{aligned} \vec\p{\bm Y} &\sim \MvNormal\p{ \vec\p{\bm X \bm\beta}, \bm\Sigma \otimes \bm I_N } \\ \Rightarrow \vec\p{\bm Y} &\sim \MvNormal\p{ (\bm I_q \otimes \bm X) \cdot \vec\p{\bm\beta}, \bm\Sigma \otimes \bm I_N } \\ \end{aligned}\]As noted by Rencher (2012) (p. 360), the BLUE for $\vec(\bm\beta)$ is
\[\begin{aligned} \vec(\hat{\bm\beta}) &= (\bm I_q \otimes (\bm X^T \bm X)^{-1} \bm X^T) \cdot \vec(\bm Y) \\ \Rightarrow \hat{\bm B} &= (\bm X^T \bm X)^{-1} \bm X^T \bm Y \end{aligned}\]with the covariance for $\vec(\hat{\bm\beta})$ being $\bm \Sigma \otimes (\bm X^T\bm X)^{-1}$. Note (happily) that the BLUE for $\bm\beta$ is not a function of $\bm \Sigma$.