Here, we only consider the codimension one case. Let $F: M\rightarrow R^{n+1}$. Nor we have the normal variation, denote like
$\phi$ is $C^\infty(R^{n+1})$, let $F(x,t)=F(x)+t\phi\nu$. Then, we have
\begin{align*}
g_{ij}(x,t)&=\langle F_i+t\phi_i\nu+t\phi\nu_i,F_j+t\phi_j\nu+t\phi\nu_j\rangle\\
&=g_{ij}(x)+t\phi\langle F_i,\nu_j\rangle+t^2\phi_i\phi_j+t^2\phi\phi_i\langle\nu,\nu_j\rangle\\
&~~+t^2\phi\langle\nu_i,F_j\rangle+t\phi\phi_j\langle \nu_i,\nu\rangle+t^2\phi^2\langle\nu_i,\nu_j\rangle
\end{align*}

Note that
\begin{align}
V(M_t)=\int_D\sqrt{det(g(x,t))}dx,
\end{align}
it follows that
\begin{align}
\frac{dV(M_t)}{dt}=\int_D\frac{1}{2\sqrt{det(g(x,t))}}\frac{\partial det(g(x,t))}{\partial t}dx
\end{align}

By elementary algebra, we have
\begin{align}
\frac{\partial det(g)}{\partial t}=det(g)trace(g^{-1}\frac{\partial g}{\partial t})
\end{align}

Note that
\begin{align*}
\frac{\partial g(x,t)}{\partial t}&=2t(\phi\phi_i\langle\nu,\nu_j\rangle+\phi\phi_j\langle \nu_i,\nu\rangle+\phi_i\phi_j+\phi^2\langle\nu_i,\nu_j\rangle)\\
&~~+\phi\langle F_i,\nu_j\rangle+\phi\langle\nu_i,F_j\rangle
\end{align*}

For $t=0$, we have
\begin{align*}
det(g(x,t))trace(g^{-1}g')\Big|_{t=0}&=\phi g^{ij}(x)\Big(\langle F_i,\nu_j\rangle+\langle\nu_i,F_j\rangle\Big)\\
&=det(g(x))\phi g^{ij}(x)\langle F_i,\nu_j\rangle++\langle\nu_i,F_j\rangle\\
&=det(g(x))\langle F_i,-h_j^kF_k\rangle+\langle-h_i^kF_k,F_j\rangle\\
&=-\phi )g^{ij}g_{ik}h_j^k-g^{ij}h_i^kg_{kj})det(g(x))\\
&=-2\phi Hdet(g(x)).
\end{align*}
Therefore,
\begin{align}
\frac{dV(M_t)}{dt}\Big|_{t=0}=-\int_{M}\phi HdV_M
\end{align}
If $\phi$ is any smooth function on $M$ with compact support, then $\frac{dV(M_t)}{dt}|_{t=0}=0$ iff $H=0$. In this case, we call $M$ is a minimal hypersurface.

Especially, if $\phi=H$, like $\frac{\partial X}{\partial t}=H\nu$, we get
\begin{align}
\frac{dV(M_t)}{dt}\Big|_{t=0}=-\int_{M}H^2dV_M\leq 0,
\end{align}
which means that under the mean curvature flow, the area of
$M_t$ is decreasing.

If we instant on assume $M_t$ is represented by $(x,u(x,t))$, by
\begin{align}
\langle\frac{dX}{dt},\nu\rangle=H=-div_{M_t}\left(\frac{(-Du,1)}{\sqrt{1+|Du|^2}}\right),
\end{align}
we have
\begin{align}
\frac{\partial u}{\partial t}=\frac{(1+|Du|^2)\Delta u-DuD^2u(Du)^T}{1+|Du|^2}.
\end{align}
This is
\begin{align}
\frac{\partial u}{\partial t}=a_{ij}(x,t)D_{ij}u,
\end{align}
where
\begin{align}
a_{ij}(x)=\delta_ij-\frac{u_iu_j}{1+|Du|^2}=g^{ij},
\end{align}
which means that the equation is strictly elliptic and the gradient estimate may play a essential in the proof of some property.

For simplicity, if we can consider the ancient similar solution in $R^n\times(-\infty,0)$, of soliton (variable separation). If $U(x,t)=\sqrt{-t}u(\frac{x}{\sqrt{-t}})$, then we have
\begin{align}
g^{ij}(Du)D_{ij}u=\frac{1}{2}x\cdot Du-\frac{1}{2}u.
\end{align}
Lu Wang's theorem assert that $U(x,t)=Du(0)\cdot x=u(x)$.

Recall that
\begin{align}
(g^{-1}(t))'=-g^{-1}g'(t)g^{-1}.
\end{align}
We will calculate $\frac{d^2V_{M_t}}{dt^2}$.

In fact, we have
\begin{align*}
\frac{\partial ^2det(g(x,t))}{\partial t^2}\Big|_{t=0}=det(g(x,t))' tr(g^{-1}g'(t))+det(g(x,t))tr((g^{-1})'g'(t)+g^{-1}g''(t))\Big|_{t=0}.
\end{align*}

Note that
\begin{align}
\partial_t\sqrt{det(g(x,t))}=\frac{\sqrt{det(g(x,t)))}}{2}tr(g^{-1}(x,t)g'(x,t)),
\end{align}
we have
\begin{align*}
\partial^2_t\sqrt{det(g(x,t))}&=\frac{\partial_t\sqrt{det(g(x,t))}}{2}tr(g^{-1}(x,t)g'(x,t))\\
&+\frac{\sqrt{det(g(x,t))}}{2}tr((g^{-1}(x,t)))'g'(x,t)+g^{-1}(x,t)g''(x,t))\\
&=\frac{1}{4\sqrt{det(g(x,t))}}\partial_t(det(g(x,t))tr\Big(g^{-1}(x,t)g'(x,t)\Big)\\
&+\frac{\sqrt{det(g(x,t))}}{2}tr\Big((g^{-1}(x,t)))'g'(x,t)+g^{-1}(x,t)g''(x,t)\Big)
\end{align*}

Firstly, we have
\begin{align*}
\partial_t(det(g(x,t))) tr(g^{-1}g'(x,t))|_{t=0}=4\phi^2 H^2det(g),
\end{align*}
that is
\begin{align*}
\frac{1}{4\sqrt{det(g(x,t))}}\partial_t(det(g(x,t))tr\Big(g^{-1}(x,t)g'(x,t)\Big)|_{t=0}=\phi^2 H^2\sqrt{det(g)}
\end{align*}
if $H\equiv0$, this term is zero.

Secondly, we have
\begin{align*}
\frac{\sqrt{det(g(x,t))}}{2}trace(g^{-1}g''(x,t))_{t=0}
=&\sqrt{det(g)}g^{ij}(\phi_i\phi_j+\phi^2\langle\nu_i,\nu_j\rangle)\\
=&(|\nabla_M\phi|^2+\phi^2|B|^2)\sqrt{det(g)},
\end{align*}
and
\begin{align*}
\frac{\sqrt{det(g(x,t))}}{2}tr((g^{-1})'g'(t))|_{t=0}&=-\sqrt{det(g)}tr(g^{-1}g'(x,0)g^-1g'(x,0))\\
&=-2\phi^2|B|^2\sqrt{det(g)}.
\end{align*}
Finally, if $H\neq 0$, we have
\begin{align}
\frac{dV_t}{dt^2}_{t=0}=\int_{M}(|\nabla _M\phi|^2+\phi^2H^2-\phi^2|B|^2)dV_M,
\end{align}
If $H\equiv0$,
we get
\begin{align}
\frac{dV_t}{dt^2}_{t=0}=\int_{M}(|\nabla _M\phi|^2-\phi^2|B|^2)dV_M,
\end{align}
we call $M$ is stable if the above quantity is nonnegative for any $\phi\in C_0^\infty(M)$.