How is local minima possible in gradient descent?
$g(x) = x^2$ is indeed a parabola and thus has just one optimum. 

However, the $\text{MSE}(\boldsymbol{x}, \boldsymbol{y}) = \sum_i (y_i - f(x_i))^2$, where $\boldsymbol{x}$ are the inputs, $\boldsymbol{y}$ the corresponding labels and the function $f$ is the model (e.g. a neural network), is not necessarily a parabola. In general, it is only a parabola if $f$ is a constant function and the sum is over one element. 

For example, suppose that $f(x_i) = c, \forall i$, where $c \in \mathbb{R}$. Then $\text{MSE}(\boldsymbol{x}, \boldsymbol{y}) = \sum_i (y_i - c)^2$ will only change as a function of one variable, $\boldsymbol{y}$, as in the case of $g(x) = x^2$, where $g$ is a function of one variable, $x$. In that case, $(y_i - c)^2$ will just be a shifted version (either to the right or left depending on the sign of $c$) of $y_i^2$, so, for simplicity, let's ignore $c$. So, in the case $f$ is a constant function, then $\text{MSE}(\boldsymbol{x}, \boldsymbol{y}) = \sum_i y_i^2$, which is a sum of parabolas $y_i^2$, which is called a paraboloid. In this case, the paraboloid corresponding to $\text{MSE}(\boldsymbol{x}, \boldsymbol{y}) = \sum_i y_i^2$ will only have one optimum, just like a parabola. Furthermore, if the sum is just over one $y_i$, that is, $\text{MSE}(\boldsymbol{x}, \boldsymbol{y}) = \sum_i y_i^2 = y^2$ (where $\boldsymbol{y} = y$), then the MSE becomes a parabola. 

In other cases, the MSE might not be a parabola or have just one optimum. For example, suppose that $f(x) = x^2$, $y_i = 1$ ($\forall i$), then $h(x) = (1 - x^2)^2$  looks as follows



which has two minima (and one maximum): at $x=0$ and $x=1$. We can find the two minima of this function $h$ using calculus: $h'(x) = -4x(1 - x^2)$, which becomes zero when $x=0$ and $x=1$. 

In this case, we only considered one term of the sum. If we considered the sum of terms of the form of $h$, then we could even have more "complicated" functions.

To conclude, given that $f$ can be arbitrarily complex, then also $\text{MSE}(\boldsymbol{x}, \boldsymbol{y})$, which is a function of $f$, can also become arbitrarily complex and have multiple minima. Given that neural networks can implement arbitrarily complex functions, then $\text{MSE}(\boldsymbol{x}, \boldsymbol{y})$ can easily have multiple minima. Moreover, the function $f$ (e.g. the neural network) changes during the training phase, which might introduce more complexity, in terms of which functions the MSE can be and thus which (and how many) optima it can have.