In the Langrangian framework without any doubt. If you are assuming an Eulerian approach the total derivative of the density has to be zero. In other words, you can have an incompresible flow with no constant fluid density as long as the local variation of the density (time derivative) cancels the convective variation ( D (dens)/Dt = - density · div v ). Take in mind that D(dens)/Dt =partial derivative of dens with respect to time + v times the gradient of density.
If the density changes in a flow are considerable (=/> 5%) that when your M=0.3 limit kicks in. For eg. for air (gamma = 1.4) This can be simply derived from isentropic relations putting the ratio of p0/p as (100/95) and solve for M.
Yes. There are many approach that totally ignore the change of density. It's a typical problem in pressure-velocity solver. SIMPLE like method for example.
One thing more I would like to focus that an incompressible flow doesn't mean that density will be uniform in flow domain. But it means that density of a fluid element should not change in time when it moves in space (path lines).
If you are solving Energy Equation along with Navier Stokes,
Simplest form will be :
dT/dt + udT/dx + udT/dy = kappa*(d2 T / dx2 +d2 T / dy2) +dissipation term (if using)