Why Normality assumption in linear regression My question is very simple: why we choose normal as the distribution that error term follows in the assumption of linear regression? Why we don't choose others like uniform, t or whatever?
 A: We do choose other error distributions. You can in many cases do so fairly easily; if you are using maximum likelihood estimation, this will change the loss function. This is certainly done in practice.
Laplace (double exponential errors) correspond to least absolute deviations regression/$L_1$ regression (which numerous posts on site discuss). Regressions with t-errors are occasionally used (in some cases because they're more robust to gross errors), though they can have a disadvantage -- the likelihood (and therefore the negative of the loss) can have multiple modes.
Uniform errors correspond to an $L_\infty$ loss (minimize the maximum deviation); such regression is sometimes called Chebyshev approximation (though beware, since there's another thing with essentially the same name). Again, this is sometimes done (indeed for simple regression and smallish data sets with bounded errors with constant spread the fit is often easy enough to find by hand, directly on a plot, though in practice you can use linear programming methods, or other algorithms; indeed, $L_\infty$ and $L_1$ regression problems are duals of each other, which can lead to sometimes convenient shortcuts for some problems).
In fact, here's an example of a "uniform error" model fitted to data by hand:

It's easy to identify (by sliding a straightedge toward the data) that the four marked points are the only candidates for being in the active set; three of them will actually form the active set (and a little checking soon identifies which three lead to the narrowest band that encompassess all the data). The line at the center of that band (marked in red) is then the maximum likelihood estimate of the line. 
Many other choices of model are possible and quite a few have been used in practice.
Note that if you have additive, independent, constant-spread errors with a density of the form $k\,\exp(-c.g(\varepsilon))$, maximizing the likelihood will correspond to minimizing $\sum_i g(e_i)$, where $e_i$ is the $i$th residual.
However, there are a variety of reasons that least squares is a popular choice, many of which don't require any assumption of normality.
A: When working with those hypothesis, squared-erros based regression and maximum likelihood provide you the same solution. You are also capable of getting simple F-tests for coefficient significance, as well as confidence intervals for your predictions.
In conclusion, the reason why we often choose normal distribution is its properties, which often make things easy. It is also not a very restrictive assumption, as many other types of data will behaive "kind-of-normally"
Anyway, as mentioned in a previous answer, there are possibilities to define regression models for other distributions. The normal just happens to be the most recurrent one
A: Glen_b has explained nicely that OLS  regression can be generalized (maximizing likelihood instead of minimizing sum of squares) and we do choose other distributions.
However, why is the normal distribution chosen so often?
The reason is that the normal distribution occurs in many places naturally. It is a bit the same like we often see the golden ratio or the Fibonacci numbers occurring "spontaneously" at various places in nature.
The normal distribution is the limiting distribution for a sum of variables with finite variance (or less strict restrictions are possible as well). And, without taking the limit, it is also a good approximation for a sum of a finite number of variables. So, because many observed errors occur as a sum of many little unobserved errors, the normal distribution is a good approximation.
See also here Importance of normal distribution
where Galton's bean machines show the principle intuitively

A: The normal/Gaussian assumption is often used because it is the most computationally convenient choice. Computing the maximum likelihood estimate of the regression coefficients is a quadratic minimization problem, which can be solved using pure linear algebra. Other choices of noise distributions yield more complicated optimization problems which typically have to be solved numerically. In particular, the problem may be non-convex, yielding additional complications. 
Normality is not necessarily a good assumption in general. The normal distribution has very light tails, and this makes the regression estimate quite sensitive to outliers. Alternatives such as the Laplace or Student's t distributions are often superior if measurement data contain outliers. 
See Peter Huber's seminal book Robust Statistics for more information. 
A: Why we don't choose other distributions?—we do.
Regression means modeling a continuous values given a set of inputs.  Consider training examples consisting of a target scalar $y_i \in \mathbb R$ and an input vector $x_i \in \mathbb R^n$.  Let the prediction of the target given $x_i$ be
$$\hat y_i = w^\intercal x_i.$$
The surprisal loss is usually the most sensible loss:
$$L = -\log P(y_i \mid x_i).$$
You can think of linear regression as using a normal density with fixed variance in the above equation:
$$L = -\log P(y_i \mid x_i) \propto (y_i - \hat y_i)^2.$$
This leads to the weight update:
$$\nabla_w L = (\hat y_i - y_i)x_i $$

In general, if you use another exponential family distribution, this model is called a generalized linear model.  The different distribution corresponds to a different density, but it can be more easily formalized by changing the prediction, the weight, and the target.
The weight is changed to a matrix $W \in \mathbb R^{n\times k}$.  The prediction is changed to
$$\hat u_i \triangleq \nabla g(W x_i)$$
where $\nabla g: \mathbb R^k \to \mathbb R^k$ is called the link function or gradient log-normalizer.  And, the target $y_i$ is changed to a vector called sufficient statistics $u_i = T(y_i) \in \mathbb R^k$.
Each link function and sufficient statistics corresponds to a different distributional assumption, which is what your question is about.  To see why, let's look at a continuous-valued exponential family's density function with natural parameters $\eta$:
$$f(z) = h(z)\exp(\eta^\intercal T(z) - g(\eta)).$$
Let the natural parameters $\eta$ be $w^\intercal x_i$, and evaluate the density at the observed target $z = y_i$.  Then, the loss gradient is
$$\begin{align}
\nabla_W L &= \nabla_W -\log f(x) \\
&= (\nabla g(W x_i)) x_i^\intercal - T(y_i) x_i^\intercal \\
&= (\hat u_i - u_i) x_i^\intercal
\end{align},$$
which has the same nice form as linear regression.

As far as I know, the gradient log-normalizer can be any monotonic, analytic function, and any monotonic, analytic function is the gradient log-normalizer of some exponential family.
