As mentioned in the comments, the question is about the expected predictive performance of a given, fitted model on new test data drawn from the data generating process (DGP). In this context, the mean squared error (MSE) is the expected squared prediction error, where the expectation is taken over the distribution of the test data, and the fitted parameters are held fixed.
According to the DGP defined in the question, let $x \sim \mathcal{N}(\vec{0}, \Sigma)$ be a new test point with corresponding response $y \sim \mathcal{N}(x \beta + \beta_0, \sigma^2)$. The predicted response is $\hat{y} = x \hat{\beta} + \hat{\beta}_0$. Under these conditions, the MSE is:
$$MSE(y, \hat{y}) = \sigma^2
+ (\beta-\hat{\beta})^T \Sigma (\beta-\hat{\beta})
+ (\beta_0-\hat{\beta}_0)^2$$
This holds regardless of how the parameters were fit (e.g. using OLS or something entirely different). The derivation is below.
MSE for a fixed test point
The MSE for a given/fixed test point $x$ is the expected squared prediction error, where the expectation is taken w.r.t. the conditional distribution of the response, given the test point:
$$MSE(y, \hat{y} \mid x) \ = \
E_{y|x} \big[ (y - \hat{y})^2 \big]$$
$$= E_{y|x} \big[ y^2 - 2 y \hat{y} + \hat{y}^2 \big]$$
By linearity of expectation:
$$= E_{y|x}[y^2] - 2 \hat{y} E_{y|x}[y] + \hat{y}^2$$
Recall that $\operatorname{Var}(y|x) = E_{y|x}[y^2] - E_{y|x}[y]^2$:
$$MSE(y, \hat{y} \mid x) \ = \
\operatorname{Var}(y|x)
+ E_{y|x}[y]^2
- 2 \hat{y} E_{y|x}[y]
+ \hat{y}^2$$
Factorize:
$$= \operatorname{Var}(y|x)
+ \Big( \hat{y} - E_{y|x}[y] \Big)^2$$
According to the DGP, the conditional mean of of the response is $E_{y|x}[y] = \beta^T x + \beta_0$ and the conditional variance is $\operatorname{Var}(y|x) = \sigma^2$. Substitute these in, along with the the predicted response $\hat{y} = \hat{\beta}^T x + \hat{\beta}_0$:
$$MSE(y, \hat{y} \mid x) \ = \
\sigma^2 + \Big(
(\hat{\beta}^T x + \hat{\beta}_0)
- (\beta^T x + \beta_0)
\Big)^2$$
Marginal MSE
The overall (marginal) MSE is obtained by averaging the MSE for a fixed test point over all possible test points, weighted by the probability density of each:
$$MSE(y, \hat{y}) = E_x \Big[ MSE(y, \hat{y} \mid x) \Big]$$
Plug in our expression for $MSE(y, \hat{y} \mid x)$ and crank through some algebra:
$$= E_x \Big[
\sigma^2
+ (\hat{\beta} - \beta)^T x x^T (\hat{\beta} - \beta)
+ 2 (\hat{\beta}_0 - \beta_0) (\hat{\beta} - \beta)^T x
+ (\hat{\beta}_0 - \beta_0)^2
\Big]$$
By linearity of expectationn:
$$= \sigma^2
+ (\beta - \hat{\beta})^T E_x[x x^T] (\beta - \hat{\beta})
+ 2 (\beta_0 - \hat{\beta}_0) (\beta - \hat{\beta})^T E_x[x]
+ (\beta_0 - \hat{\beta}_0)^2$$
According to the DGP, the distribution over test points has mean zero: $E_x[x] = \vec{0}$. Thus, $E_x[x x^T] = \Sigma$ is the covariance matrix of the predictors and the marginal MSE is:
$$MSE(y, \hat{y}) \ = \ \sigma^2
+ (\beta - \hat{\beta})^T \Sigma (\beta - \hat{\beta})
+ (\beta_0 - \hat{\beta}_0)^2$$