This is not automatic (yet), but it's not difficult to work around it by simply adding another rule for how $V$ should vary. Something like
\partial{#}::PartialDerivative;
\delta{#}::Accent;
g^{\mu\nu}::Symmetric;
S:= \int{ \sqrt{-g} ( 1/2 g^{\mu\nu}\partial_{\mu}{\phi} \partial_{\nu}{\phi} + V) }{x};
rl:= \phi -> \delta{\phi}, V -> V' \delta{\phi};
vary(S, rl)
sort_product(S)
canonicalise(_)
will produce
$$\int \sqrt{-g}(\partial_{\mu} \delta{\phi} \partial_{\nu} \phi g^{\mu\nu} + V^\prime \delta{\phi}) {\rm d}x$$
after which you can do
distribute(_)
integrate_by_parts(_, $\delta{\phi}$)
factor_out(_, $\delta{\phi}$);
to get things in an equation-of-motion form.