The quartic
\begin{equation}
f = x^4 - 7x^3y + 16x^2y^2 - 7xy^3 + y^4 + 2x^3 - 8x^2y - 8xy^2 + 2y^3 + x^2 - xy + y^2
\end{equation}
has isolated local minima at the five places
\begin{equation}
(0, 0), (0, -1), (-1, 0), (2, 5), (5, 2).
\end{equation}
In order to prove the claim, one checks that the partial derivatives $\frac{\partial f}{\partial x}$ and $\frac{\partial f}{\partial y}$ vanish in these points, and that the Hessian $\begin{pmatrix}\frac{\partial^2 f}{\partial x^2}&\frac{\partial^2 f}{\partial x\partial y}\\\frac{\partial^2 f}{\partial x\partial y}&\frac{\partial^2 f}{\partial y^2}\end{pmatrix}$ is positive definite in these five points. A sufficient (and necessary) condition for that is that the traces and the determinants are positive.

If one carries out the calculation by hand, then using the symmetry $f(x,y)=f(y,x)$ allows us to check only $3$ points. Alternatively, the following Sage code verifies the example:

    R.<x, y> = QQ[]
    f = (x^4 - 7*x^3*y + 16*x^2*y^2 - 7*x*y^3 + y^4 + 2*x^3 - 
         8*x^2*y - 8*x*y^2 + 2*y^3 + x^2 - x*y + y^2)
    fx, fy = f.derivative(x), f.derivative(y)
    critical_points = ideal([fx, fy]).variety()
    hesse = matrix(2, 2, [fx.derivative(x), fx.derivative(y),
                          fy.derivative(x), fy.derivative(y)])
    minima = []
    for cp in critical_points:
        h = hesse.substitute(cp)
        if h.trace() > 0 and h.det() > 0:
            minima.append(cp)
    print(len(minima) == 5)

*Remark 1:* As pointed out by Jap88 in the comments, this example is difficult to visualize. The reason is that three saddle points are very close to three minima. The red and blue lines are the curves given by $\frac{\partial f}{\partial x}=0$ and $\frac{\partial f}{\partial y}=0$, respectively. The second and third picture zoom into those difficult places:[![enter image description here][1]][1]

*Remark 2:* If we measure the size of a polynomial by the sum of the absolute values of its coefficients, and look for quartic polynomials $f(x,y)$ with integral coefficients, $5$ integral local minima one of which is $(0,0)$, and which is symmetric in $x$ and $y$, then I have numerical evidence that the example in the answer is optimal.

*Remark 3:* Conjecture 5.4 in [this paper by Durfee et.al][2] predicts that a quartic has at most $4$ isolated local minima. So this is a counterexample.

*Remark 4*: Theorem 3.1.6. of the paper [Critical points of real polynomials ...][3] by Shustin seems to claim that $5$ minima are possible. However, the "proof" given in Section 3.13 merely says "Assume that $d \le 4$. Then all the index distributions except [...] can be easily realized as it was explained above". I don't know what the author means by "above", as all the previous cases handle degrees $d\ge5$. 


  [1]: https://i.sstatic.net/x6O5j.png
  [2]: http://www.jstor.org/stable/2324459
  [3]: https://www.researchgate.net/profile/Eugenii-Shustin-2/publication/267655438_Critical_points_of_real_polynomials_subdivisions_of_Newton_polyhedra_and_topology_of_real_algebraic_hypersurfaces/links/593ec9480f7e9bf167c5c86c/Critical-points-of-real-polynomials-subdivisions-of-Newton-polyhedra-and-topology-of-real-algebraic-hypersurfaces.pdf?origin=publicationDetail&_sg%5B0%5D=q9Oh1_m3jOK6GbcVqAfUCbUZQcUhzBdCSuiy3C4-M4tptZBRj9TUFvwMlxungIYZGBJCN448_Ym-sIGtnLa7Tg.4tpQfI3YawFUZmb-f-azz2eQ4vhkYXKkq72yY3dxP10EjQoQVfzrjiDaqocMD5tauPGPA_EHk65Awgyzqi5H6A&_sg%5B1%5D=SGk8BZqoH6nDDCqLNnHvxuM-t9dfzYMqDyjU0p1TCYgtjQarC-uAAejb8-TZRRJc3tDQhpv1GlfrHxvmiaWwho-UDx5azx9mjQ7Oto0W3pHr.4tpQfI3YawFUZmb-f-azz2eQ4vhkYXKkq72yY3dxP10EjQoQVfzrjiDaqocMD5tauPGPA_EHk65Awgyzqi5H6A&_iepl=&_rtd=eyJjb250ZW50SW50ZW50IjoibWFpbkl0ZW0ifQ%3D%3D