Why are the Least-Squares and Maximum-Likelihood methods of regression not equivalent when the errors are not normally distributed? Title says it all. I understand that the Least-Squares and Maximum-Likelihood will give the same result for regression coefficients if the model's errors are normally distributed. But, what happens if the errors are not normally distributed? Why are the two methods no longer equivalent?
 A: Because the MLE is derived from the assumption of residual normally distributed.
Note that 
$$
\text{min}_\beta~~ \|X \beta - y \|^2
$$
Has no probabilistic meaning: just find the $\beta$ that minimize the squared loss function. Everything is deterministic, and no random components in there.
Where the concept of probability and likelihood come, is we assume 
$$
y=X\beta + \epsilon
$$
Where we are considering $y$ as a random variable, and $\epsilon$ is normally distributed.
A: The least squares and the maximum (gaussian) likelihood fit are always equivalent. That is, they are minimized by the same set of coefficients.
Changing the assumption on the errors changes your likelihood function (maximizing the likelihood of a model is equivalent to maximizing the likelihood of the error term), and hence the function will no longer be minimized by the same set of coefficients.
So in practice the two are the same, but in theory, when you maximize a different likelihood, you will get to a different answer than Least-squares
A: Short Answer
The probability density of a multivariate Gaussian distributed variable $x=(x_1, x_2,...,x_n)$, with mean $\mu=(\mu_1,\mu_2,...,\mu_n)$ is related to the square of the euclidean distance between the mean and the variable ($\vert \mu-x \vert_2^2$), or in other words the sum of squares.

Long Answer
If you multiply multiple Gaussian distributions for your $n$ errors, where you assume equal deviations, then you get a sum of squares.
$$ \begin{array}
\mathcal{L(\mu_j,x_{ij})} = P(x_{ij} \vert \mu_j) & =\prod_{i=1}^n \frac{1}{\sqrt{2 \pi \sigma^2}} exp\left[-\frac{(x_{ij}-\mu_i)^2}{2\sigma^2}\right] \\
&= \left(\frac{1}{\sqrt{2 \pi \sigma^2}} \right)^n exp \left[  -\frac{\sum_{i=1}^n(x_{ij}-\mu_i)^2}{2\sigma^2}\right]
\end{array}$$
or in the convenient logarithmic form:
$$ 
\log\left(\mathcal{L(\mu_j,x_{ij})} \right) = n \log \left( \frac{1}{\sqrt{2 \pi \sigma^2}} \right) -\frac{1}{2\sigma^2} \sum_{i=1}^n(x_{ij}-\mu_j)^2
$$
So optimizing the $\mu$ to minimize the sum of squares is equal to maximizing the (log) likelihood (ie. the product of multiple Gaussian distributions, or the multivariate Gaussian distribution).
It is this nested square of the difference $(\mu-x)$ inside exponential structure, $exp\left[ (x_i-\mu)^2 \right]$, which other distributions do not have.

Compare for instance with the case for Poisson distributions
$$\log(\mathcal{L}) = \log \left( \prod\frac{\mu_j^{x_{ij}}}{x_{ij}!} exp \left[ -\mu_j \right] \right) = -\sum \mu_j -\sum log(x_{ij}!) + \sum log(\mu_j) x_{ij} $$
which has a maximum when the following is minimized:
$$\sum \mu_j -log(\mu_j) x_{ij}$$
which is a different beast.

In addition (history)
The history of the normal distribution (ignoring deMoivre getting to this distribution as an approximation for the binomial distribution) is actually as the discovery of the distribution that makes the MLE correspond to the least squares method (rather than the the least squares method being a method that can express the MLE of the normal distribution, first came the least squares method, second came the Gaussian distribution)
Note that Gauss, connecting the 'method of maximum likelihood' with the 'the method of least squares', came up with the 'Gaussian distribution', $e^{-x^2}$ , as the sole distribution of errors that leads us to make this connection between the two methods.
From Charles Henry Davis' translation (Theory of the motion of the heavenly bodies moving about the sun in conic sections. A translation of Gauss's "Theoria motus," with an appendix) ...
Gauss defines:

Accordingly, the probability to be assigned to each error $\Delta$ will be expressed by a function of $\Delta$ which we shall denote by $\psi \Delta$.
(Italization done by me)

And continues (in section 177 pp. 258):

... whence it is readily inferred that $\frac{\psi^\prime\Delta}{\Delta}$ must be a constant quantity. which we will denote by $k$. Hence we have $$\text{log } \psi \Delta = \frac{1}{2} k \Delta \Delta + \text{Constant}$$ $$\psi \Delta = x e^{\frac{1}{2}k \Delta \Delta}$$ denoting the base of the hyperbolic logarithms by $e$ and assuming $$\text{Constant} = \log x$$

ending up (after normalization and realizing $k<0$) in

$$\psi \Delta = \frac{h}{\sqrt{\pi}} e^{-hh\Delta \Delta}
$$

A: A concrete example: Suppose we take a simple error function p(1)=.9, p(-9) =.10 . If we take two points, then LS is just going to take the line through them. ML, on the other hand, is going to assume that both points are one unit too high, and thus will take the line through the points shifted down on unit.
