First of all, gauge invariance is a property of a
classical theory. We expect that it is preserved under quantization, i.e. that it still holds in quantum field theory.
If this is valid then every mathematical object, which is a scalar w.r.t the action of the gauge group classically, is translated into a scalar in quantum field theory. So in essence the classical algebraic property holds in quantum field theory.
In principle, it is possible that this is
not valid, i.e. that quantization
breaks gauge invariance. That means that quantization introduced a so-called
gauge anomaly. So proving gauge invariance means proving the absence of gauge-anomalies.
I do not know any example where the gauge invariance is broken by an anomaly
and where the resulting quantum field theory is still consistent. Therefore, all consistent quantum gauge theories must be gauge invariant and must preserve gauge invariance of mathematical objects which are gauge invariant classically i.e. algebraically.
How to prove gauge invariance is nicely demonstrated in the case of the so-called Adler-Bell-Jackiw- or chiral anomaly
https://en.wikipedia.org/wiki/Chiral_anomaly
Here gauge invariance can be preserved but that implies that global chiral symmetry has to become anomalous. You may look for the original papers or have a look in some textbooks like Ryder (as far as I remember).
In the standard model (to be precise, in the electro-weak sector) the chiral symmetry is essentially the gauge symmetry, due to the fact that different chiralities of fermions couple differently to the gauge fields. That means that gauge symmetry is anomalous per fermion species, but that summing over all fermion species (involved in specific processes i.e. loops) the individual contributions add up to zero. Therefore, the overall gauge symmetry is preserved.
The proof for one specific correlator simply means to calculate this correlator summing over all contributions from different fields. The proof the “the whole theory” means that you have to have a scheme how a finite (hopefully small :-) subset of correlators in low order of perturbation theory ensures this property for all correlators up to all orders. This is closely related with techniques known from proving (perturbative) renormalizability. A modern approach is provided by Fujikawa.
The question whether gauge invariant correlators are always constructed from gauge invariant operators is related to the question whether all formally (algebraically) gauge invariant correlators can always be constructed from gauge invariant operators acting on the gauge invariant vacuum. If you study something like
\sum_{a,b}\langle \phi_a|\mathcal{O}_{ab} |\psi_b\rangle
that would mean that the equation
\sum_{a,b}\langle \phi_a|\mathcal{O}_{ab} |\psi_b\rangle = \langle 0| \sum_{a,b} \mathcal{X}_a \,\mathcal{O}_{ab} \, \mathcal{Y}_b |0\rangle
must hold.Hope that helps.