Comment: the following is a somewhat convoluted way of deriving the
Euler-Lagrange equation using Clairaut's theorem for the volume functional and
some standard, albeit not simpler, variation formulas (all is $C^{\infty}$ in
the following). Let $g_{0}$ be a Riemannian metric and let $v$ be a symmetric
$2$-tensor. Let $g_{0,s}$ satisfy $\frac{\partial}{\partial s}|_{s=0}
g_{0,s}=v$ and $g_{0,0}=g_{0}$. With $g_{0,s}$ as initial data, let $g_{t,s}$,
$t\in\lbrack0,\varepsilon)$, solve the Ricci-Hamilton-DeTurck flow
$\frac{\partial}{\partial t}g_{t,s}=-2\operatorname{Ric}{}_{g_{t,s}
}+\mathcal{L}_{W_{t,s}}g_{t,s}$, where $W_{t,s}=\operatorname{tr}^{1,2}
{}_{g_{t,s}}(\nabla_{g_{t,s}}-\nabla_{g_{t,0}})$. Note that $\frac{\partial
}{\partial t}g_{t,0}=-2\operatorname{Ric}{}_{g_{t,0}}$. Let $v_{t,s}
=\frac{\partial}{\partial s}g_{t,s}$. We have $\frac{\partial^{2}}{\partial
s\partial t}|_{s=0}g_{t,s}=\Delta_{L}v_{t,0}$, which equals $\frac
{\partial^{2}}{\partial t\partial s}|_{s=0}g_{t,s}=\frac{\partial}{\partial
t}v_{t,0}$ (Lichnerowicz Laplacian heat equation). We compute $\frac{\partial
}{\partial t}\operatorname{Vol}(g_{t,s})=\frac{1}{2}\int\operatorname{tr}
{}_{g_{t,s}}(\frac{\partial}{\partial t}g_{t,s})d\mu_{g_{t,s}}=-\int
R_{g_{t,s}}d\mu_{g_{t,s}}$ since $\int\operatorname{tr}{}_{g_{t,s}
}(\mathcal{L}_{W_{t,s}}g_{t,s})d\mu_{g_{t,s}}=0$ by the divergence theorem.
Now
\begin{align*}
\frac{\partial}{\partial s}|_{s=0}\int R_{g_{t,s}}d\mu_{g_{t,s}} &
=-\frac{\partial^{2}}{\partial s\partial t}|_{s=0}\operatorname{Vol}
(g_{t,s})=-\frac{\partial^{2}}{\partial t\partial s}|_{s=0}\operatorname{Vol}
(g_{t,s})\\
& =-\frac{1}{2}\frac{\partial}{\partial t}\int\operatorname{tr}{}_{g_{t,0}
}(v_{t,0})d\mu_{g_{t,0}}\\
& =\int\langle-\operatorname{Ric}{}_{g_{t,0}}+\frac{g_{t,0}}{2}R_{g_{t,0}
},v_{t,0}\rangle_{g_{t,0}}d\mu_{g_{g_{t,0}}}
\end{align*}
since $\int\operatorname{tr}{}_{g_{t,0}}(\frac{\partial}{\partial t}
v_{t,0})d\mu_{g_{t,0}}=\int\operatorname{tr}{}_{g_{t,0}}(\Delta_{L}
v_{t,0})d\mu_{g_{t,0}}=\int\Delta_{g_{t,0}}(\operatorname{tr}{}_{g_{t,0}
}(v_{t,0}))d\mu_{g_{t,0}}=0$. Finally, take $t=0$.
December 18, 2013. The notion of volume, in various guises, occurs throughout
the study of Ricci flow, especially in Perelman's work. Now, per unit increase
in scale $t$, the volume form of a metric changes with velocity $\frac
{\partial}{\partial t}d\mu=-Rd\mu$. By Clairaut's theorem, the variation of
$-Rd\mu$ is equal to the change per unit increase in scale of the variation of
the volume form, i.e.,
$$
\frac{\partial}{\partial t}(\frac{\operatorname{tr}_{g}v}{2}d\mu_{g})=\left(
\langle\operatorname{Ric}-\frac{1}{2}Rg,v\rangle+\operatorname{div}
(\frac{\nabla\operatorname{tr}v}{2})\right) d\mu.
$$
In the $f$-warped or entropy version of this, we have $\frac{\partial
}{\partial t}(fe^{-f}d\mu)=(-R-\Delta f)e^{-f}d\mu$ under $\frac{\partial
}{\partial t}g=-2(\operatorname{Ric}+\nabla^{2}f)$ and $\frac{\partial
f}{\partial t}=-R-\Delta f$. Integrating this yields that $\mathcal{N}
\doteqdot\int_{\mathcal{M}}fe^{-f}d\mu$ satisfies $-\frac{d\mathcal{N}}
{dt}=\mathcal{F}\doteqdot\int(R+|\nabla f|^{2})e^{-f}d\mu$. If $\frac
{\partial}{\partial s}g=v$ and $\frac{\partial f}{\partial s}=\frac
{\operatorname{tr}_{g}v}{2}$, then the variation of the energy integrand is
\begin{align*}
& \frac{\partial}{\partial s}((-R-\Delta f)e^{-f}d\mu)\\
& =\left( (-L(v,\nabla f)+2\langle\operatorname{Ric}+\nabla^{2}
f,v\rangle)e^{-f}+\operatorname{div}(e^{-f}\{\frac{\nabla\operatorname{tr}
v}{2}-v(\nabla f)\})\right) d\mu\doteqdot A,
\end{align*}
where $L(v,X)=\operatorname{div}^{2}v+\langle\operatorname{Ric},v\rangle
-2\langle\operatorname{div}v,X\rangle+v(X,X)$ is the linear trace Harnack
quadratic. On the other hand, $\frac{\partial}{\partial s}(fe^{-f}d\mu
)=\frac{\operatorname{tr}_{g}v}{2}e^{-f}d\mu$. So Perelman's version is
$\frac{\partial}{\partial t}(\frac{\operatorname{tr}_{g}v}{2}e^{-f}d\mu
)=\frac{\partial^{2}}{\partial s\partial t}(fe^{-f}d\mu)=A$, using Clairaut's
theorem. Note that integration by parts gives $\int L(v,\nabla f)e^{-f}
d\mu=\int\langle\operatorname{Ric}+\nabla^{2}f,v\rangle e^{-f}d\mu$, from
which one obtains Perelman's energy variation formula.
In Section 6.2 of arXiv:0211159 Perelman argues that the $\mathcal{W}$-entropy
(i.e., $\mathcal{F}$ with scaling) integrand is a warped scalar curvature. So,
without scaling (i.e., $\tau$), we would have the correspondences
$\mathcal{N}\sim\operatorname{Vol}$ and $\mathcal{F}\sim\int Rd\mu$, which is
also clear from taking $f=\operatorname{const}\neq0$ as a special case.
However, in 6.2, Perelman's volume is essentially $\int e^{-f}d\mu$, which is
constant under the above variations.
Best Answer
This can be found in Besse "Einstein Manifolds", in chapter 4.
The idea is to use Koszul formula for the Levi-Civitta connection to compute the derivative of the curvature with respect to the metric. Bianchi identities also help.