Since the OLS model assumptions are that $y_i=\beta^\top x_i + \epsilon_i$ with $\epsilon_i\stackrel{\text{i.i.d.}}{\sim}\mathcal N(0,\sigma^2)$, $\text{MS}_{\text{residual}}:=\frac 1 n \sum_i(y_i-\hat{\beta}^\top x_i)^2$ is (only) an approximation of the variance $\sigma^2$ for two reasons
1. $\hat \beta\approx\beta\in\mathbb{R}^d$ (holds in expectation)
The OLS estimator $\hat\beta$ is unbiased as $\mathbb{E}\left[\hat \beta\right]=\beta$ and the expected error converges to 0 in the limit since $\mathbb{E}\left[||\hat\beta-\beta||^2\right]=\frac {\sigma^2d}{n}$
2. $\frac 1 n \sum_i(y_i-\beta^\top x_i)^2=\frac 1 n \sum_i\epsilon_i^2 \approx\sigma^2$ (holds in limit)
By the law of large numbers, as $n \rightarrow \infty$ we have $\frac 1 n\sum_i\epsilon_i^2\rightarrow\mathbb{E}\left[\epsilon_i^2\right]=\sigma^2$
Search "Statistical guarantees for ordinary least squares" for proofs and such