Consider the standard regression model with i.i.d. observations $(X_i,Y_i)$ for $i=1,2,\dots,n$:
$$ Y_i = \beta_0 + \beta_1 X_{i} + \varepsilon_i, \quad \quad i = 1,2,\dots,n, $$ where the regressors $X_i$ are considered to be random variables as opposed to fixed observations, and the errors are normally distributed conditional on the regressors and have fixed variance.
Suppose we solve this model using ordinary least squares and obtain estimated coefficients $\hat \beta_0, \hat \beta_1$.
Now define $n$ new random variables $$ Z_i = \hat \beta_0 + \hat \beta_1 X_{i}, \quad \quad i=1,2,\dots,n. $$
How do we calculate $\text{Var}[\sum_{i=1}^n Z_i]$? I'm not sure if $Z_i$ are independent because they are constructed using $\hat \beta_0$ and $\hat \beta_1$ which makes it seem like the $Z_i$ could be dependent on each other?
Note: I want to treat $\hat \beta_0$ and $\hat \beta_1$ as random. Wikipedia says these estimates are normally distributed since the errors are normally distributed.