How does the distribution of the error term affect the distribution of the response? So when I assume that the error terms are normally distributed in a linear regression, what does it mean for the response variable, $y$?
 A: We invent the error term by imposing a fictitious model on real data; the distribution of the error term does not affect the distribution of the response.
We often assume that the error is distributed normally and thus try to construct the model such that our estimated residuals are normally distributed. This can be difficult for some distributions of $y$. In these cases, I suppose you could say that the distribution of the response affects the error term.
A: Maybe I'm off but I think we ought to be wondering about $f(y|\beta, X)$, which is how I read the OP. In the very simplest case of linear regression if your model is $y=X\beta + \epsilon$ then the only stochastic component in your model is the error term. As such it determines the sampling distribution of $y$. If $\epsilon\sim N(0, \sigma^2I)$ then $y|X, \beta\sim N(X\beta, \sigma^2I)$. What @Aniko says is certainly true of $f(y)$ (marginally over $X, \beta$), however. So as it stands the question is slightly vague.
A: If you write out the response as
$$\bf{y}=m+e$$
Where $\bf{m}$ is the "model" (the prediction for $\bf{y}$) and $\bf{e}$ is the "errors", then this can be re-arranged to indicate $\bf{y}-m=e$.  So assigning a distribution for the errors is the same thing as indicating the ways your model is incomplete.  To put it another way is that it indicates to what extent you don't know why the observed response was the value that it actually was, and not what the model predicted.  If you knew your model was perfect, then you would assign a probability distribution with all of its mass on zero for the errors.  Assigning a $N(0,\sigma^{2})$ basically says that the errors are small in units of $\sigma$.  The idea is that the model predictions tend to be "wrong" by similar amounts for different observations, and is "about right" on the scale of $\sigma$.  As a contrast, an alternative assignment is $Cauchy(0,\gamma)$ which says that most of the errors are small, but some errors are quite large - the model has the occasional "blunder" or "shocker" in terms of predicting the response.
In a sense the error distribution is more closely linked to the model than to the response.  This can be seen from the non-identifiability of the above equation, for if both $\bf{m}$ and $\bf{e}$ are unknown then adding an arbitrary vector to $\bf{m}$ and subtracting it from $\bf{e}$ leads to the same value of $\bf{y}$, $\bf{y}=m+e=(m+b)+(e-b)=m'+e'$.  The assignment of an error distribution and a model equation basically says which arbitrary vectors are more plausible than others.
A: The short answer is that you cannot conclude anything about the distribution of $y$, because it depends on the distribution of the $x$'s and the strength and shape of the relationship. More formally, $y$ will have a "mixture of normals" distribution, which in practice can be pretty much anything.
Here are two extreme examples to illustrate this:


*

*Suppose there are only two possible $x$ values, 0 an 1, and $y = 10x + N(0,1)$. Then $y$ will have a strongly bimodal distribution with bumps at 0 and 10.

*Now assume the same relationship, but let $x$ be uniformly distributed on the 0-1 interval with lots of values. Then $y$ will be almost uniformly distributed over the 0-10 interval (with some half-normal tails at the edges).


In fact, since every distribution can be approximated arbitrarily well with mixture of normals, you can really get any distribution for $y$.
