I tried to convince myself it is Diff(M), but I failed. Most books say Bianchi Identities reduce the independent equations in Einstein's equations by 4, therefore there are some redundancies in the metric variables. As a result, there could be many solutions that correspond to one physical state. There are two arguments which show the gauge group is Diff(M). Hole argument: a diffeomorphism that preserves a Cauchy surface and moves other points will cause the Cauchy data to not have a unique evolution, therefore two metrics related by a diffeomorphism must correspond to one physical state. Objection: the set of all diffeomorphisms that preserve a Cauchy surface is not Diff(M), there exists diffeomorphisms that move all points around. Dirac constraint analysis: GR has two first class constraints, one corresponds to 3D diffeomorphism of the hypersurface (Diffeomorphism constraint), one corresponds to diffeomorphism that moves forward in time (Hamiltonian constraint). Objection: the two constraints generate diffeomorphisms on a hypersurface (at an instant of time), not the whole of M. So I believe the gauge group of GR is not Diff(M), it should be a subset of it.