Your argument only applies to finite systems (otherwise the energy is ill-defined) and there are no phase transitions in finite systems. So, there is no contradiction there.
Moreover, your argument only applies when both $h=0$ (no magnetic field) and you use free or periodic boundary conditions. Indeed, were it not the case, then you would not have symmetry under spin-flip.
Now, consider a system in the box $\{-n,\dots,n\}^d$ with, say, $+$ boundary condition (that is, all spins in the exterior boundary of the box are fixed to $+1$). Let us denote the corresponding probability measure by $\mu_{n,\beta}^+$ and the associated expectation by $\langle\cdot\rangle_{n,\beta}^+$.
Then (assuming that $d\geq 2$), one can (rather easily) show, using for instance Peierls' argument, that at low enough temperatures, the expected value of the central spin $\sigma_0$ is positive: there exist $\epsilon>0$ and $\beta_0$ (both independent of $n$) such that, for all $\beta>\beta_0$,
$$
\langle\sigma_0\rangle_{n,\beta}^+ \geq \epsilon.
$$
In the same way, one shows that, for all $\beta>\beta_0$,
$$
\langle\sigma_0\rangle_{n,\beta}^- \leq - \epsilon,
$$
for a system with $-$ boundary condition.
We now want to define probability measures on the set of all infinite configurations (that is, configurations of all spins in $\mathbb{Z}^d)$. I will not enter into too much detail here. One way to do that is to take the thermodynamic limit. That is, we would like to define a measure $\mu^+_\beta$ as the limit of $\mu^+_{n,\beta}$ as $n\to\infty$. The precise sense in which this limit is taken is the following one: for any local observable $f$ (that is, any observable depending only on the values taken by finitely many spins), we want convergence of the expectation of $f$:
$$
\langle f \rangle_\beta^+ = \lim_{n\to\infty} \langle f \rangle_{n,\beta}^+.
$$
One can show, using correlation inequalities, that the limit indeed exists in this sense. Moreover, in view of the above, for all $\beta>\beta_0$,
$$
\langle \sigma_0 \rangle_\beta^+ = \lim_{n\to\infty} \langle \sigma_0 \rangle_{n,\beta}^+ \geq \epsilon.
$$
One can do the same starting with the $-$ boundary condition and define a measure $\mu^-_\beta$ as the limit of the measures $\mu^-_{n,\beta}$ and we'll have, for all $\beta>\beta_0$,
$$
\langle \sigma_0 \rangle_\beta^- \leq -\epsilon.
$$
In particular, the two measures $\mu^+_\beta$ and $\mu^-_\beta$ cannot coincide (since the expectation of $\sigma_0$ is different under these two measures!). You have thus shown that your system can exist in two different phases when there is no magnetic field and the temperature is low enough. In the phase described by $\mu^+_\beta$, the magnetization is positive, while it is negative in the phase described by $\mu^-_\beta$.
Of course, you might also have considered the limit of measures with free (or periodic) boundary conditions $\mu^\varnothing_\beta$ and have concluded that, for all $\beta$,
$$
\langle \sigma_0\rangle_\beta^\varnothing = 0.
$$
However, the measure $\mu^\varnothing_\beta$ does not describe a pure phase. In fact,
$$
\mu_\beta^\varnothing = \frac12\mu^+_\beta + \frac12\mu^-_\beta .
$$
Pure phases are important for several reasons. First, these are the only ones in which macroscopic observables take on deterministic values. Second, they contain all the interesting physics, since any other Gibbs measure $\mu$ can be written as a convex combination of pure phases (as we did above for $\mu_\beta^\varnothing$). In particular, if you sample a configuration with $\mu$, then you'll obtain a configuration that is typical of one of the pure phases (with a probability corresponding to the corresponding coefficient in the convex decomposition; for instance, using $\mu_\beta^\varnothing$, you'd obtain a configuration typical of $\mu^+_\beta$ with probability $1/2$).
(Pure phases possess additional remarkable properties, but this would take us too far, so I'll only discuss this if requested explicitly.)
Let me briefly describe an alternative way to proceed. Rather than introducing boundary conditions that break the symmetry, you can continue to work with, say, periodic boundary condition, but introduce a magnetic field $h$. Denote the corresponding measure $\mu_{n,\beta,h}^{\rm per}$.
Then, one can again take the limit as $n\to\infty$ and obtain a limiting measure $\mu_{\beta,h}$. This measure can be shown to be unique as long as $h\neq 0$, in the sense that the limit does not depend on the boundary condition used. Moreover, one has that
$$
\lim_{h\downarrow 0} \mu_{\beta,h} = \mu^+_\beta
$$
and
$$
\lim_{h\uparrow 0} \mu_{\beta,h} = \mu^-_\beta.
$$
So, the two measures obtained before, describing the pure phases of the (infinite-volume) Ising model, correspond precisely to the phases you obtain by setting a positive (resp. negative) magnetic field and decreasing (resp. increasing) it to $0$.
Combined with the discussion above, this explains how the magnetization can have a discontinuity at $h=0$ at low temperatures.
To conclude (finally!), let me just mention that it is possible to construct infinite-volume Gibbs measures (such as the measures $\mu_\beta^+$ and $\mu^-_\beta$ described above) directly in infinite-volume, without taking limits of finite-volume measures. This is interesting because this avoids any explicit symmetry breaking! I discussed this in another answer.
