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
∫√−g(∂μδϕ∂νϕgμν+V′δϕ)dx
after which you can do
distribute(_)
integrate_by_parts(_, $\delta{\phi}$)
factor_out(_, $\delta{\phi}$);
to get things in an equation-of-motion form.