First, we need the expression for the propagator. Let's start with the expression for the Greens function
\begin{align}
\square G = \delta^4(x)\,.
\end{align}
We now perform a Fourier transform, so the right hand side of the equation becomes
\begin{align}
\int{\rm d}^4 x e^{{\rm i}(\omega t - \vec k\cdot\vec x)}\delta^4(x) = 1\,.
\end{align}
The left hand side of the equation is
\begin{align}
\square G(x) = (\partial_t^2 - \nabla^2)G(x)\to (-\omega^2 + (\vec k)^2)G(k)\,.
\end{align}
Thus,
\begin{align}
G(k) = \frac{-1}{\omega^2 - (\vec k)^2}\,.
\end{align}
The retarded propagator is obtained by the following prescription, $\omega\to \omega + {\rm i}\epsilon$. Thus, we can perform the inverse Fourier transform to obtain the position space representation of $G$:
\begin{align}
G_\text{ret}(r) = \int\frac{{\rm d}^4 k}{(2\pi)^4} \frac{e^{-{\rm i}(\omega t - \vec k\cdot\vec x)}}{(\omega + {\rm i}\epsilon)^2 - (\vec k)^2}\,.
\end{align}
Importantly, the poles of this function are in the lower half-plane. Therefore, the integral must close in the lower half plane in order to be non-zero (see the residue theorem for an explanation). When $t > 0$, the coefficient of $\omega$ is negative, so we close the integral in the lower half plane. Similarly, when $t < 0$, we close the integral in the upper half plane. The latter case does not contribute, which we express with a Heaviside function
\begin{align}
G_\text{ret} &= \theta(t)2\pi{\rm i}\int\frac{{\rm d}^3 k}{(2\pi)^4}\left( \frac{e^{-{\rm i}(-k t - \vec k\cdot\vec x)}}{2k} + \frac{e^{-{\rm i}(k t - \vec k\cdot\vec x)}}{2k}\right)\,,\\
&=\theta(t)\int\frac{{\rm d}^3 k}{(2\pi)^3}e^{{\rm i} k x}\frac{\sin kt}{k}
\end{align}
We now take the angular integral
\begin{align}
G_\text{ret} &= \theta(t)\int \frac{4\pi k^2{\rm d} k}{(2\pi)^3}j_0(k r)\frac{\sin kt}{k}\,,\\
&=\frac{\theta(t)}{r}\int \frac{4\pi {\rm d} k}{(2\pi)^3}\sin(k r)\sin(k t)\,.
\end{align}
We now observe that this is just an expression of the $\delta$-function
\begin{align}
G_\text{ret} &=\frac{\theta(t)}{4\pi r}\delta(t - r)\,.
\end{align}
Now that we have the correct expression for the Greens function, let's apply the D'Alembertian to it and see that this is indeed the Greens function.
\begin{align}
\square\frac{\theta(t)}{4\pi r}\delta(t - r) &= \left(\partial_t^2 - \nabla^2\right)\frac{\theta(t)}{4\pi r}\delta(t - r)\,,\\
&=\frac{1}{4\pi r}\partial_t^2\theta(t)\delta(t - r) - \theta(t)\nabla^2\frac{1}{4\pi r}\delta(t - r)\,,\\
&=\frac{1}{4\pi r}\partial_t(\delta(t - r)\delta(t) + \theta(t)\delta'(t - r)) - \theta(t)\left(\delta(t - r)\nabla^2\frac{1}{4\pi r} + \frac{1}{4\pi r}\nabla^2\delta(t - r) + 2\partial_r\delta(t - r)\partial_r\frac{1}{4\pi r}\right)\,,\\
&=\frac{1}{4\pi r}(2\delta'(t - r)\delta(t)+\delta(t - r)\delta'(t) +\theta(t)\delta''(t - r)) - \theta(t)\left(-\delta(t - r)\delta(r) + \frac{1}{4\pi r}\delta''(t - r)- \frac{2}{4\pi r^2}\delta'(t - r) + 2\delta'(t - r)\frac{1}{4\pi r^2}\right)\,,\\
&=\frac{1}{4\pi r}(2\delta'(t - r)\delta(t)+\delta(t - r)\delta'(t) +\theta(t)\delta''(t - r)) - \theta(t)\left(-\delta(t - r)\delta(r) + \frac{1}{4\pi r}\delta''(t - r)\right)\,,\\
&=\frac{1}{4\pi r}(2\delta'(t - r)\delta(t)+\delta(t - r)\delta'(t)) - \theta(t)\left(-\delta(t - r)\delta(r)\right)\,,\\
&=\theta(t)\delta(t - r)\delta(r)\,.
\end{align}
Note, the identity that $\nabla^2 1/r = -4\pi\delta(r)$ is well explained elsewhere (I'll provide a link later).
