The answer from Tim Mak is perfect and brings some further readings about what happens out of the usual distributional assumptions. But I would like to bring another view or intuition on the matter.
You are telling me that our coefficient estimates will be the same in both situations, as long as you make the assumption that the data is normal, even if it is not normal?
Yes and no matter what you assume really. You can always get those parameters, using the very same method of least-square or MLE, with arbitrarily false assumptions. At worst your model fit will be bad (and wrong). People often fit a linear model and look afterwards at the error distribution: if it is not gaussian you might find that a better model is needed.
Intuitively
First, let's remember that the MLE is a bayesian account of linear regression and treats the data themselves as random variables. It gives a probabilistic model of the data, and hence we ought to start with a distribution somewhere (the assumption on the error is a good place to start).
On the other hand, a least-square estimation is a simple approach for parameter estimation, we only assume a model and a measure of errors, parameters are the variable. Gauss himself, by trying to incorporate a probability density on the error of such model came up with the "gaussian" as we can read on the Wikipedia page:
He then turned the problem around by asking what form the density should have and what method of estimation should be used to get the arithmetic mean as estimate of the location parameter. In this attempt, he invented the normal distribution.
In one case (LSE) we find the parameter $\beta$ such as the euclidian distance (equivalently the mean squared error) is minimised. Here we decide to measure the error as the square distance.
In the other (MLE) we expect an error distributed normally around the line $y = \beta x$, which means an error that occurs more or less likely depending on the squared distance between data pairs. Then the parameter estimate that maximizes the probability of our data under this "error expectation" (or assumption) are the same parameters that we found by deciding to take the square distance as our error metric.
It simply is due to how you measure errors or discrepancy between your model (here we say a linear one such as $\beta x$) and your observables ($y$). Both LSE and MLE rely on measuring this as related to the square distance between those, although from a different starting point or way of thinking.
It can be helpful to remember that gaussianity in general "pairs well" with 2-norm errors (or also that laplacian distribution "pairs" with 1-norm errors...).
More maths
If might feel natural to come back to the equations anyway...
The key element is simply that both the classical Least Square Estimate and the Maximum Likelihood Estimation for linear regression (with gaussianity assumption for errors) rely on the square distance between each data point pairs $(y_i, \beta x_i)$ as an error metric or measure of discrepancy between both datasets.
- For the least-square perspective, it feels natural to use the euclidian distance to measure how far apart our model $\beta x_i$ and the observable $y_i$ are. So our cost function will be of the form:
$$
\sum_i(y_i - \beta x_i)^2
$$
- For the MLE, it is our most uninformative choice to assume a normally distributed error of data around the model mean for a given standard deviation (which is like saying that the error should remain evenly distributed around 0, a nice euclidian symmetry as offered by the square distance we see in the exponential $e^{-1/2(\frac{\epsilon}{\sigma})^2}$). So all we insist on is having a model of $y$ such as:
$$
y = \beta x + \epsilon \\
\text{with}:~~\epsilon \sim \mathcal{N}(0, \sigma)
$$
And this gives us
$$
\implies ~~ y \sim \mathcal{N} (\beta x, \sigma) = \frac{1}{\sqrt{2\pi\sigma}} . e^{-\frac{1}{2}\left(\frac{(y-\beta x)}{\sigma}\right)^2}
$$
Meaning that once again, if we were to find the most likely $\beta$ for such i.i.d. data we end up minimizing $$\sum_i(y_i - \beta x_i)^2$$
