Is the linear combination of least square estimators in linear regression normally distributed? I am reading Mathematical Statistics with Applications by Wackerly et al. (7th edition). In Chapter 11, the book discusses linear models and least squares, specifically $Y = \beta_0+\beta_1x+\epsilon$ where $E(\epsilon)=0$ and $Var(\epsilon)=\sigma^2$, with $E(Y)=\beta_0+\beta_1x$ being deterministic. 
The book derives $\hat\beta_0=\bar Y - \hat\beta_1\bar x$ and $\hat\beta_1=\frac{\sum(x_i-\bar x)(Y_i-\bar Y)}{\sum(x_i - \bar x)^2}$ as estimators of $\beta_0$ and $\beta_1$ respectively from minimizing the SSE. The book states if $\epsilon$ is normally distributed, then each $Y_i$ from an independent sample is normally distributed and  $\hat\beta_0$ and $\hat\beta_1$ are also normally distributed since they are linear combinations of $Y_i$. I am assuming this is true based on this proof? 
However, the book later goes further to state that $\hat\beta_0$ and $\hat\beta_1$ being normal implies that a linear combination $\hat\theta = a_0\hat\beta_0 + a_1\hat\beta_1$ is also normal. This confuses me since $Cov(\hat\beta_0 ,\hat\beta_1)=\frac{-\bar x\sigma^2}{\sum(x_i-\bar x)^2}$ may not be zero, which means for $\hat\theta$ to be normal, the joint distribution of $\hat\beta_0$ and $\hat\beta_1$ must be normal (per this proof) but the book makes no reference to the joint distribution. Is it always true that $\hat\theta$ is normally distributed? If so, how is it shown?
 A: The answer is yes. This is because least squares estimates themselves are jointly normal distributed which can be seen by applying the linear transformation theorem for the multivariate normal distribution (Proof):
$$
x \sim \mathcal{N}(\mu,\Sigma) \quad \Rightarrow \quad y = Ax + b \sim \mathcal{N}(A\mu + b, A \Sigma A^\mathrm{T})
$$
Using matrix notation for linear regression, the distribution of the data is:
$$
y = X\beta + \varepsilon, \; \varepsilon \sim \mathcal{N}(0, \sigma^2 I_n)
$$
$$
\Rightarrow \quad y \sim \mathcal{N}(X\beta, \sigma^2 I_n)
$$
If we apply ordinary least squares (Proof), the distribution of the estimates is:
$$
\hat{\beta} = (X^\mathrm{T}X)^{-1} X^\mathrm{T}y
$$
$$
\Rightarrow \hat{\beta} \sim \mathcal{N}\left((X^\mathrm{T}X)^{-1} X^\mathrm{T}X\beta, \sigma^2 (X^\mathrm{T}X)^{-1} X^\mathrm{T} X (X^\mathrm{T}X)^{-1} \right) = \mathcal{N}\left( \beta, \sigma^2 (X^\mathrm{T}X)^{-1} \right)
$$
Then, a linear combination of the least squares estimates is distributed as:
$$
\Rightarrow c^\mathrm{T} \hat{\beta} \sim \mathcal{N}\left( c^\mathrm{T} \beta, \sigma^2 c^\mathrm{T} (X^\mathrm{T}X)^{-1} c \right)
$$
In your special case (a.k.a. "simple linear regression"), the quantities would be:
$$
X = \begin{bmatrix} 1_n & x \end{bmatrix}, \; \beta = \begin{bmatrix} \beta_0 \\ \beta_1 \end{bmatrix}, \; c = \begin{bmatrix} a_0 \\ a_1 \end{bmatrix}, \; c^\mathrm{T} \hat{\beta} = \hat{\theta}
$$
$$
( \text{and of course} \quad y = Y, \; \varepsilon = \epsilon )
$$
