The connection is made through a specific procedure: namely, ordinary least squares. When you suppose a response $Y$ is a random variable related to other variables $X$ in the form
$$Y = X\beta + \varepsilon$$
where $\varepsilon$ is assumed to have a joint Normal distribution (plus some other assumptions that don't matter here), then the least squares solution is
$$\hat\beta = (X^\prime X)^{-}X^\prime Y = (X^\prime X)^{-}X^\prime (X\beta + \epsilon) = \beta + J\varepsilon$$
for a matrix $J = (X^\prime X)^{-}X^\prime $ that is computed only from $X.$ The residuals, or estimated errors, are the differences
$$\hat\varepsilon = Y - \hat Y = (X\beta + \varepsilon) - X\hat\beta = \varepsilon - XJ\varepsilon = (\mathbb I - H)\varepsilon$$
where $H = XJ = X(X^\prime X)^{-}X^\prime$ is the "hat matrix." This exhibits the residuals as linear combinations of $\varepsilon$ (with coefficients given by $\mathbb I - H$). Linear combinations of jointly Normal variables are Normal, QED.