Is it really proper to say Ward identity is a consequence of gauge invariance?

This answer partially disagrees with Motl's. The crucial point is to consider the difference between the abelian and non-abelian case. I totally agree with Motl's answer in the non-abelian event — where these identities are usually denominated Slavnov-Taylor's rather than Ward's, so that I will refer to the abelian case.

First, a few words about terminology: Ward identities are the quantum counterpart to (first and second) Noether's theorem in classical physics. They apply to both global and gauge symmetries. However, the term is often reserved for the $U(1)$ gauge symmetry in QED. In the case of gauge symmetries, Ward identities yield real identities, such as $k^{mu}mathcal M_{mu}=0$, where $mathcal M_{mu}$ is defined by $mathcal M=epsilon_{mu},mathcal M^{mu}$, in QED, that tell us that photon's polarizations parallel to photon's propagation don't contribute to scattering amplitudes. In the case of global symmetries, however, Ward identities reflect properties of the theory. For example, the S-matrix of a Lorentz invariant theory is also Lorentz invariant or the number of particles minus antiparticles in the initial state is the same as in the final state in a theory with a global (independent of the point in space-time) $U(1)$ phase invariance.

Let's study the case of a massive vectorial field minimally coupled to a conserved current:

$$mathcal L=-{1over 4},F^2+{a^2over 2}A^2+i,barPsidisplaystyle{not}D, Psi - {m^2over 2}barPsiPsi =-{1over 4},F^2+{a^2over 2}A^2+i,barPsidisplaystyle{not}partial , Psi - {m^2over 2}barPsiPsi-e,A_{mu},j^{mu}$$

Note that this theory has a global phase invariance $Psirightarrow e^{-itheta},Psi$, with a Noether current

$$j^{mu}={barPsi, gamma^{mu}},Psi$$

such that (classically) $partial_{mu},j^{mu}=0$. Apart from this symmetry, it is well-known that the Lagrangian above is equivalent to a theory: i)that doesn't have an explicit mass term for the vectorial field. ii) that contains a scalar field (a Higgs-like field) with a different from zero vacuum expectation value, which spontaneously break a $U(1)$ gauge symmetry (this symmetry is not the gauged $U(1)$ global symmetry mentioned previously). The equivalence is in the limit where vacuum expectation value goes to infinity and the coupling between the vectorial field and the Higgs-like scalar goes to zero. Since one has to take this last limit, the charge cannot be quantized and therefore the $U(1)$ gauge symmetry must be topologically equivalent to the addition of real numbers rather than the multiplication of complex numbers with unit modulus (a circumference). The difference between both groups is only topological (does this mean then that the difference is irrelevant in the following?). This mechanism is due to Stueckelberg and I will summarize it at the end of this answer.

In a process in which there is a massive vectorial particle in the initial or final state, the LSZ reductio formula gives:

$$langle i,|,f ranglesim epsilon _{mu}int d^4x,e^{-ikcdot x}, left(eta^{munu}(partial ^2-a^2)-partial^{mu}partial^{nu}right)...langle 0|mathcal{T}A_{nu}(x)...|0rangle$$

From the Lagrangian above, the following classical equations of motion may be obtained

$$left(eta^{munu}(partial ^2-a^2)-partial^{mu}partial^{nu}right)A_{nu}=ej^{mu}$$

Then, quantumly,

$$left(eta^{munu}(partial ^2-a^2)-partial^{mu}partial^{nu}right)langle 0|mathcal{T}A_{nu}(x)...|0rangle = e,langle 0|mathcal{T}j^{mu}(x)...|0rangle + text{contact terms, which don't contribute to the S-matrix}$$

And therefore

$$langle i,|,f ranglesim epsilon _{mu}int d^4x,e^{-ikcdot x},...langle 0|mathcal{T}j^{mu}(x)...|0rangle +text{contact terms, which don't contribute}sim epsilon_{mu}mathcal{M}^{mu}$$

If one replaces $epsilon_{mu}$ with $k_{mu}$, one obtains

$$k_{mu}mathcal{M}^{mu}sim k _{mu}int d^4x,e^{-ikcdot x},...langle 0|mathcal{T}j^{mu}(x)...|0rangle$$

Making use of $k_{mu}sim partial_{mu},,e^{-ikcdot x}$, integrating by parts, and getting ride of the surface term (the plane wave is an idealization, what one actually has is a wave packet that goes to zero in the spatial infinity), one gets

$$k_{mu}mathcal{M}^{mu}sim int d^4x,e^{-ikcdot x},..., partial_{mu},langle 0|mathcal{T}j^{mu}(x)...|0rangle$$

One can now use the Ward identity for the global $Psirightarrow e^{-itheta},Psi$ symmetry (classically $partial_{mu},j^{mu}=0$ over solutions of the matter, $Psi$, equations of motion)

$$partial_{mu}, langle 0|mathcal{T}j^{mu}(x)...|0rangle = text{contact terms, which don't contribute to the S-matrix}$$

And hence

$$k^{mu}mathcal M_{mu}=0$$

same as in the massless case.

Note that in this derivation, it has been crucial that the explicit mass term for the vectorial field doesn't break the global $U(1)$ symmetry. This is also related to the fact that the explicit mass term for the vectorial field can be obtained through a Higgs-like mechanism connected with a hidden (the Higgs-like field decouples from the rest of the theory) $U(1)$ gauge symmetry.

A more careful calculation should include counterterms in the interacting theory, however I think that this is the same as in the massless case. We can think of the fields and parameters in this answer as bare fields and parameters.

Stueckelberg mechanism

Consider the following Lagrangian

$$mathcal L=-{1over 4},F^2+|dphi|^2+mu^2,|phi|^2-lambda, (phi^*phi)^2$$

where $d=partial - ig, B$ and $F$ is the field strength (Faraday tensor) for $B$. This Lagrangian is invariant under the gauge transformation

$$Brightarrow B + (1/g)partial alpha (x)$$ $$phirightarrow e^{ialpha(x)}phi$$

Let's take a polar parametrization for the scalar field $phi$: $phiequiv {1over sqrt{2}}rho,e^{ichi}$, thus

$$mathcal L=-{1over 4},F^2+{1over 2}rho^2,(partial_{mu}chi-g,B_{mu})^2+{1over 2}(partial rho)^2+{mu^2over 2},rho ^2- {lambdaover 4}rho^4$$

We may now make the following field redefinition $Aequiv B - (1/g)partial chi$ and noting that $F_{munu}=partial_{mu}B_{nu}-partial_{nu}B_{mu}=partial_{mu}A_{nu}-partial_{nu}A_{mu}$ is also the field strength for $A$

$$mathcal L=-{1over 4},F^2+{g^2over 2}rho^2,A^2+{1over 2}(partial rho)^2+{mu^2over 2},rho ^2-{lambdaover 4}, rho^4$$

If $rho$ has a vacuum expectation value different from zero $langle 0|rho |0rangle = v=sqrt{mu^2over lambda}$, it is then convenient to write $rho (x)=v+omega (x)$. Thus

$$mathcal L=-{1over 4},F^2+{a^2over 2},A^2+g^2,v,omega,A^2+{g^2over 2},omega ^2,A^2+{1over 2}(partial omega)^2-{mu^2over 2},omega ^2-lambda,vomega^3-{lambdaover 4}, omega^4+{v^4,lambda^2over 4}$$

where $aequiv gtimes v$. If we now take the limit $grightarrow 0$, $vrightarrow infty$, keeping the product, $a$, constant, we get

$$mathcal L=-{1over 4},F^2+{a^2over 2},A^2+{1over 2}(partial omega)^2-{mu^2over 2},omega ^2-lambda,vomega^3-{lambdaover 4}, omega^4+{v^4,lambda^2over 4}$$

that is, all the interactions terms between $A$ and $omega$ disappear so that $omega$ becomes an auto-interacting field with infinite mass that is decoupled from the rest of the theory, and therefore it doesn't play any role. Thus, we recover the massive vectorial field with which we started.

$$mathcal L=-{1over 4},F^2+{a^2over 2},A^2$$

Note that in a non-abelian gauge theory must be non-linear terms such as $sim g A^2,partial A;$, $sim g^2 A^4$, which prevent us from taking the limit $grightarrow 0$.