We know from the exact solutions for 2D Ising model on square lattice the long range order appears bellow critical temperature, but how does this agree with the Mermin-Wagner theorem, from which we know that there is no long range order in 2D space above $T=0$?
There are no phase transitions in 1D. In 2D, the Mermin-Wagner theorem "states that continuous symmetries cannot be spontaneously broken at finite temperature in systems with sufficiently short-range interactions" as Wikipedia put it. The important thing is that Z2 Ising has a discrete symmetry, so the theorem does not apply. The other condition, that the interaction is short range is indeed fullfilled.