Variation of the Green function with respect to the metric

Consider a (closed) Riemann surface and let $G(x,y)$ be the Green function of the Laplace-Beltrami operator. We can informally identify $G$ with the two-point correlation function for the Gaussian random field:

$$G(x,y)=\left<\phi(x)\phi(y)\right>=\frac{1}{Z} \int \mathcal{D}\phi\;\phi(x)\phi(y) \exp\left(-\frac{1}{2}\int \left| \nabla\phi(z)\right|^2 dV_g(z)\right).$$ This is easy to vary with respect to the metric $g_{\mu\nu}$:

$$\delta G(x,y)=-\frac{1}{2}\int dV_g(z)\;\delta g^{ij}(z)\left<\phi(x)\phi(y) \nabla_i\phi(z)\nabla_j\phi(z)\right>_c,$$ where $\left<\right>_c$ is the connected correlation function (the disconnected diagram $x\leftrightarrow y,z\leftrightarrow z$ is cancelled by the variation of the partition function). We can use Wick's theorem to compute this, getting

$$\boxed{\frac{\delta G(x,y)}{\delta g^{ij}(z)}=-\nabla_{[i} G(x,z)\nabla_{j]} G(z,y)}$$

(derivatives w.r.t. $z$). This formula looks wickedly similar to the Hadamard variation formula for the variation of the boundary of the domain in flat space. Yet I haven't been able to find any mentions of this in the mathematical literature.

Furthermore, if we define the regularized Green's function (a.k.a. the Robin function) by $$G^R(x)=\lim_{y\to x} \left(G(x,y)-\frac{1}{2\pi}\ln d(x,y)\right),$$ where $d(x,y)$ is the local geodesic distance, then I'm conjecturing the following variational formula: $$\frac{\delta G^R(x)}{\delta g^{ij}(z)}=-\nabla_{i} G(x,z)\nabla_{j} G(x,z)-\frac{1}{4\pi}\nabla_i \nabla_j G(x,z).$$ The second term is motivated by the well-known formula for conformal variations (where it is $\frac{1}{4\pi}\delta_x(z)$) and seems to be necessary to cancel the second order pole in this variation. Edit: this guess turned out to be wrong, see my answer below.

$$\frac{\delta G(x,y)}{\delta g^{\mu\nu}(z)}=(-\nabla_{[\mu}G(z,x)\nabla_{\nu]}G(z,y)+\frac{1}{2}g_{\mu\nu}(z)\nabla_\rho G(z,x)\nabla^\rho G(z,y))+\frac{1}{2V}(G(x,z)+G(z,y))g_{\mu\nu}(z).$$
The traceless part in the brackets represents the quasiconformal variation (which in complex conformal coordinates reduces to just $8\partial_z G(z,x)\partial_z G(z,y)$, as found in the literature). The last "trace" term is the conformal variation, which is not probed by the path integral because it only shows up at finite volume.
I believe this formula to be correct based on the independence of conformal and quasiconformal variations. The resulting tensor is also divergence-free away from $x$ and $y$, as required by general covariance. The more difficult part is proving the following conjecture for the Robin function:
$$\frac{\delta G^{R}\left(x\right)}{\delta g^{\mu\nu}\left(z\right)}=-\nabla_{\mu}G\left(z,x\right)\nabla_{\nu}G^R\left(z\right)+\frac{1}{2}g_{\mu\nu}(z)\nabla_\rho G(z,x)\nabla^\rho G^R(z)+\frac{1}{V}G\left(z,x\right)g_{\mu\nu}\left(z\right)-\frac{1}{4\pi}\nabla_{\mu}\nabla_{\nu}G\left(z_{i},z\right)-\frac{\delta(x,z)}{4\pi}g_{\mu\nu}(z)+\frac{1}{8\pi V}g_{\mu\nu}(z)+\frac{1}{2V}f_{\mu\nu},$$ where $f_{\mu\nu}$ is an unknown symmetric traceless tensor defined by $$\nabla^\mu f_{\mu\nu}=\nabla_\nu G^R,$$ as required by the divergence-free condition on the metric variation of $G^R$. This result is based on the exact formula $$G^R(z)=\frac{1}{4\pi}\int G\left(x,y\right)R\left(y\right)\mathrm{d}V(y)+\frac{1}{V}\zeta^{R}\left(1\right)-c,$$ where $\zeta^{R}(s)=\zeta(s)-\frac{V}{4\pi}\frac{1}{s-1}$ and $\zeta$ is the spectral zeta function of the positive Laplace-Beltrami operator, $c=\frac{1}{2\pi}\left(\gamma-\ln2\right)$ (reference below).
Steiner, Jean, A geometrical mass and its extremal properties for metrics on $S^2$, Duke Math. J. 129, No. 1, 63-86 (2005). ZBL1144.53055.