Varying a scalar field Lagrangian density

I was varying a scalar field density and I look at this term $${\cal L}~=~-\frac{1}{2}\partial _\mu\phi\partial^\mu\phi.$$

The result that I need to come is $$-\frac{1}{2}\delta(\partial _\mu\phi\partial^\mu\phi)=\nabla_\mu\nabla^\mu\phi\delta\phi.$$ But I can't find a way to do it. How the partial derivates transform into the nabla?

Reference:

1. https://arxiv.org/abs/1809.09274. This is the paper where I see it, the action is in (1) and the variation result is in (4).

Hint: The square root $$\sqrt{|g|}$$ in the measure is here important. The infinitesimal variation of the action term is
$$\delta \left( \int \! d^4x~\sqrt{|g|} \frac{1}{2}\partial_{\mu}\phi ~g^{\mu\nu}~\partial_{\nu}\phi\right) ~=~ \int \! d^4x~\sqrt{|g|} \partial_{\mu}\delta\phi ~g^{\mu\nu}~\partial_{\nu} \phi$$ $$~\stackrel{\text{int. by parts}}{\sim}~- \int \! d^4x~\delta\phi~\partial_{\mu}(\sqrt{|g|} g^{\mu\nu}\partial_{\nu} \phi) ~=:~- \int \! d^4x~\sqrt{|g|}~\delta\phi~\nabla^2\phi,$$ because of the formula $$\nabla^2 ~:=~g^{\nu\mu}\nabla_\nu\nabla_\mu~=~\frac{1}{\sqrt{|g|}}\partial_{\mu}\sqrt{|g|} g^{\mu\nu}\partial_{\nu}$$ for the Laplace-Beltrami operator, see e.g. my Phys.SE answer here.
It is the through the 4d-version of the 3d vector identity $$\nabla\cdot (\delta \phi \nabla \phi)=(\nabla \delta \phi) \cdot (\nabla\phi)+\delta \phi \nabla^2 \phi$$ together with an integration by parts that gives $$\delta\left\{\frac 12 \int (\nabla\phi)\cdot (\nabla\phi) d^3x\right\}$$ $$= \int (\nabla\delta \phi)\cdot (\nabla\phi) d^3x$$ $$= - \int \delta \phi \nabla^2 \phi \,d^3x$$ There is no difference really between the 3d and 4 cases except for the range of summation over the repeated indices.