homework and exercises - Schwinger-Dyson equations: Deriving an equation of motion for $langle phi(x) phi(y) rangle$

I'm considering the Euclidean Klein-Gordon theory, with action

$$S_{0}[\phi] = \frac{1}{2}\int~\mathrm d^{4}x ~\phi(x) \left[ - \partial_{x}^{2} + m^{2} \right] \phi(x).\tag{1}$$ My generating function is then given by:

$$\mathcal{Z}_{0}[J] = \int \mathcal{D}[\phi]\ \exp\left( - S_{0}[\phi] + \int ~\mathrm d^{4}x \ \phi(x) J(x)\right).\tag{2}$$

I'm supposed to derive the following equation of motion, involving the two-point function:

$$\left[ - \partial_{y}^{2} + m^{2} \right] \langle \phi(y) \phi(x) \rangle = \delta^{(4)}(y-x).\tag{3}$$

I've been told to do this by starting with the following definition of the one-point function:

$$\langle \phi(x) \rangle = \frac{1}{\mathcal{Z}_{0}[0]} \frac{\delta \mathcal{Z}_{0}[J]}{\delta J(x)} \bigg|_{J~=~0} \tag{4}$$

I'm supposed to use the invariance of the functional integration under field re-defintions; IE. if we replace $\phi(x)$ with $\phi^{\prime} = \phi(x) + \epsilon(x)$.

$\ $

I've been looking online and the usual approach that I've seen involves looking at the equality

$$\int \mathcal{D}[\phi] \exp( - S_{0}[\phi] ) \phi(x) = \int \mathcal{D}[\phi^{\prime}] \exp( - S_{0}[\phi^{\prime}] ) \phi^{\prime}(x) \tag{5}$$ and performing an expansion in $\epsilon$.

How would you do this starting with $\langle \phi(x)\rangle$?