why not using sample variance (instead of MSE) to estimate the error variance in linear regression? Assuming the true equation for Y is linear as below:
$$Y_i =\beta_1X_i +\beta_0 + \epsilon_i$$
Assuming X is fixed, then the variance of each Y is:
$$var(Y_i )=var(\epsilon_i)=\sigma^2$$
In order to estimate $\sigma ^2$, we usually use the mse:
$$\hat{\sigma}^2=\frac{1}{n-p-1}\sum_i(Y_i -\hat{Y_i})^2 $$
My question is, why can't we estimate $\sigma^2$ using the sample variance below?
$$\hat{\sigma}^2=\frac{1}{n-1}\sum_i (Y_i-\bar{Y})^2$$
The sample variance is also an unbiased estimator, so instead of predicted value we just use the sample average of Y?

$\hat\sigma=\frac{1}{N-1}\sum_{i=1}^N\left(Y_i-\bar Y\right)^2$ has nothing to do with the model, unless your model is to predict the mean $\bar Y$ every single time, no matter what the features are.
A measure of prediction quality that is unrelated to the predictions does not make much sense.
In fact, the reason for doing regression is to tighten up your estimates and do a better job of predicting the conditional mean than you would do by predicting the pooled/marginal mean every time.
Your $\hat\sigma$ does come up in $R^2$, however, as $R^2=1-\dfrac{
\sum_{i=1}^N\left(
Y_i -\hat Y_i
\right)^2
}{
(N-1)\frac{1}{N-1}\sum_{i=1}^N\left(Y_i-\bar Y\right)^2
}=1-\dfrac{
\sum_{i=1}^N\left(
Y_i -\hat Y_i
\right)^2
}{
(N-1)\hat\sigma
}$.
EDIT
Here is a quick simulation to show that the two equations give dramatically different results, even for a "large" sample size of a million.
set.seed(2022)
N <- 1000000
x <- runif(N, 0, 12)
y <- x + rnorm(N)
L <- lm(y ~ x)
(1 / (N - 2)) * sum((y - predict(L))^2) # I get 0.9992724.
(1 / (N - 1)) * sum((y - mean(y))^2) # I get 12.98154.

EDIT 2
"Sigma" has two meanings. You are right that $\sigma^2$ can be used to denote the constant variance of the errors. This would related to the conditional distribution, which is our object of interest in regression, and $\frac{1}{n-p-1}\sum_i(Y_i -\hat{Y_i})^2 $ is unbiased for this conditional variance.
However, $\sigma^2$ also could apply just to $Y$, marginally/pooled, in which case, $\frac{1}{n-1}\sum_i(Y_i -\bar{Y})^2 $ is an unbiased estimator for this unconditional/marginal/pooled variance.
If you want to differentiate between them, perhaps consider using subscripts like $\sigma^2_{Y\vert X}$ for the conditional variance and $\sigma^2_{Y}$ for the unconditional variance, but they refer to different random variables that have different variances.
