How is Y Normally Distributed in Linear Regression I understand that $y_{i}$s are normally distributed because we assume that the residual is normally distributed which seems a decent assumption. 
Question: Does that mean $Y$ is also normally distributed? In other words, if I plot a histogram of $Y$ values will they have a unimodal Gaussian shape? 
Confusion: We generally apply linear regression on problems that involve predicting a continuous variable. Given the wide applicability of linear regression, I am finding it little difficult to believe that $Y$ tend to be normally distributed for all of them. 
I believe that in case of Logistic Regression, I can say that $Y$ will be Bernoulli distributed, and $y_{i}$s will have the same distribution. In other words, I can say that $Y$ has a Bernoulli distribution with $P(y_{i}=1) = \theta$
 A: I will try to kill your confusion in the one dimensional case. 
The proof with $k$ regressors is simply an extension and does not affect the proof below.
We have ONLY the following hypotesis:


*

*$Y = \beta X + e $

*$e \sim N(0,  \sigma^2)$

*$X$ and $e$ are independent
Given these hypotesis you can prove that $Y\mid X \sim N(\beta x, \sigma^2)$ 
BUT you CANNOT prove that $Y$ is Normal. In order to do that you need one additional assumption  about the JOINT DISTRIBUTION of $X$ and $e$ that is $X$ and $e$ are JOINTLY Normal distributed. This of course implies that $X$ and $Y$ are JOINTLY Normal distributed as well.
We are not making any assumption about the distribution of $X$, and I want to point that $X$ is RANDOM unlike many people said that $X$ is NON RANDOM. That is not true and you'll do yourself no favor if you think that.
Consequently $X$ has a probability distribution. That being said let's begin with the proof:
$f_{Y \mid X}(y \mid x) = \dfrac{f_{YX}(y,x)}{f_{X}(x)}$
$f_{YX}(y,x)=f_{YX}(\beta X + e,x) = f_{eX}(y- \beta X,x)$
the last equality is true because it is a simple bivariate transformation 
The independence of $e$ and $X$ implies that 
$f_{eX}(y- \beta X,x) = f_e(y- \beta X)f_X(x)$
Finally, since $e \sim N(0,  \sigma^2)$
$f_{Y \mid X}(y \mid x) = f_e(y- \beta X) = \dfrac{1}{2 \sigma \sqrt{(\pi)}}exp{\bigg(\dfrac{1}{2} \Big(\dfrac{y -\beta x}{\sigma} \Big)^2 \bigg)} $
This conclude the proof since the last equality implies that.
$ Y \mid X \sim N(\beta x, \sigma^2)$
To summarize: given the hypotesis in the OP, if you don't know the joint distribution of $X$ and $Y$ you CANNOT draw any conclusion about the UNCONDITIONAL of $Y$ BUT you can determine the CONDITIONAL DISTRIBUTION of $Y \mid X$ independently of the MARGINAL of $X$
A: $$y_i = a + \sum b_{j} x_{ij} + \epsilon_i$$
You can have a multitude of distribution types for $y_i$, $x_{ij}$ and $\epsilon_i$ in this linear regression. While the three, $y_i$, $x_{ij}$, and $\epsilon_i$, relate to each other, I guess that it is clear that when we only restrict/assume a single one of them, the $\epsilon_i \sim \mathcal{N(0,\sigma^2)}$ (which is the special case of linear regression called OLS, ordinary least squares, regression), that the $y_i$ and $x_i$ can still have some independent degree of variation between each other and do not need to be normal distributed like the error term $\epsilon$.

Example: Y can be modeled as a quadratic function of a normal distributed variable X plus some normal distributed error term. In this case you can clearly fit Y using OLS (fitting to a quadratic function) but Y is not normal distributed (since it is related to a sum of variables in which one of them is a quadratic function of a normal distribution.
an example image found in another stackexchange post:
https://stats.stackexchange.com/a/30981/164061
 

The variations/generalizations, of OLS, that you encounter in regression (such as GLM) assume 


*

*different (more general) distributions of the error term $\epsilon$ 

*as well apply some non linear function of the ensemble of the linear
part $\sum b_j x_{ij}$ (note this is different from the $x_{ij}$ which might be non-linear functions of some measured set of parameters such as in quadratic fits, these cases do not take  non-linear manipulation of the entire linear sum) 


The presence of these options might be confusing and generate the idea that OLS is for normal distributed $Y$ and GLM and those other cases for non-normal distributed $Y$. But, the case is more nuanced and OLS can also apply to non-normal distributed $Y$ (and more specifically is just about normal distributed error-term).
A: Regarding the statement marked question, for linear regression, if the vector of residuals $e \sim N(0, \sigma^2 I)$ then since $y = X \beta + e$ and $X \beta$ is non-random $y \sim N(X \beta, \sigma^2 I)$ so, yes, $y$ is normally distributed, as well.
Regarding the histogram, normally histograms are made from values that are independent and identically distributed but the observed $y_i$ values under the linear model have different distributions as their means differ so that does not fulfill the requirement for a histogram (except for the trivial linear model consisting of intercept only).
Regarding the statement marked confusion, you might want to look up the Gauss Markov theorem which provides an optimality condition for least squares without the assumption of normality.
Regarding logistic regression, that is not a special case of linear regression.
Logistic regression is an example of a generalized linear model (GLM) and GLMs form a family that is a superset which include both the linear model and logistic regression.
Linear regression (referred to in the subject of the post and above in this answer) refers to regression with a normally distributed response variable.  The predictor variables and coefficients are fixed (i.e. non-random) and the residuals are normally distributed as well.   In R one uses the lm function to analyze such models.
Generalized linear regression (GLM) is a superset of linear regression.  The assumptions are somewhat similar to linear models but now the dependent variable belongs to an exponential family (not to be confused with the exponential distribution).  That family includes the normal, binomial, exponential, poisson and other distributions.  The mean of the response variable is related via a link function to a linear function of the independent variables (which is where the linearity comes in).  If the dependent variable is normally distributed and the link function is the identity function then GLMs reduce to linear models.  If the binomial distribution is assumed then the model is referred to as binomial logistic regression.  This can encompass 0/1 data (binary logistic regression).  In R one uses the glm function to analyze GLMs.  Note that the general linear model is not the same as the generalized linear model (GLM) so be careful when reading accounts of this. 
There are other directions in which one can generalize but normally these are not what one is referring to when one refers to linear regression and these have separate names.  For example, one could assume that the $\beta$ values are random (random effects model) or that some are fixed and some are random (mixed effects models).  Such models are fit by the lme4 package in R.   These are closely related to Bayesian models.  Other models include nonlinear regression models where the residuals are still normal and independent but $X \beta$ is replaced with a nonlinear function of $\beta$.  The nls function in R fits such models.  Measurement error models ($X$ is random), generalized mixed models (GLMs where coefficients are partly or entirely random), time series models (notably auto-correlation and ARIMA dependence structures among the residuals), multivariate response models ($Y$ is a matrix) and other models may have a form similar to linear regression models as well.  
Regarding references, McCullagh and Nelder, Generalized Linear Models (second edition) is useful and has separate chapters for different types of data.  For example, Chapter 4 is about binary data.
Another possibility is to find a book with a chapter on GLMs.  Here are some at a variety of levels. Chapter 2 of N. Wood, Generalized Additive Models, is about GLMs and discusses both theory as well as examples using R code -- the book as a whole is about GAMs but this chapter is about GLMs.  Zuur et al, Mixed Effect Models and Extensions in Ecology with R is a relatively non-mathematical book that discusses GLMs in Chapters 8, 9 and 10 together with R code.  McCulloch et al, Generalized, Linear and Mixed Models provides a terse discussion of GLMs in Chapter 5.
Note
Note that a number of those commenting have claimed that $X$ is random. If you do make that assumption then of course the answer would change; however, as you were asking How is Y Normally Distributed in Linear Regression not only is that not the standard assumption but that route will lead to it not being distributed normal whereas the standard assumptions do lead to it being normal so they fulfill the statement of the question as summarized in its subject.
Another interpretation (not mine) is that you are asking about the distribution of the observed values $y_i$ as represented by their histogram. In that case as pointed out in comments under the question it has already been answered here 
What if residuals are normally distributed, but y is not? and the question is a duplicate.  Again, be aware that that will not lead to the requested explanation of how Y is normal.
EDIT:  Have expanded the discussion at the end and added a paragraph on the histogram.
Have moved comments to Note above as comments were getting too long and contentious.
A: Actually, your y scores need not be normally distributed. Even when y departs from normality, it is entirely possible for the residuals to assume normality, which is one of the assumptions of linear regression. 
The distinction with logistic regression (or poisson or other), however, is a somewhat different issue. There you are specifying the mapping (or link) between changes in X and y scores. 
