The "one phrase answer" is "divergence theorem".


Slightly wordier but a bit formally (for ease of typing I write $dz = dx~dv$ for the volume on phase space)

$$ \partial_t \int f^p dz = \int \partial_t f^p dz $$

Next, 

$$ 0 = \int \nabla_x \cdot (vf^p) dz $$ 

assuming $f$ decays suitably fast at infinity, and similarly

$$ 0 = \int \nabla_v \cdot (Ff^p) dz $$

Now, 

$$ \partial_t f^p + \nabla_x \cdot (vf^p) + \nabla_v \cdot (F f^p) = p f^{p-1} \left[ \partial_t + v\cdot \nabla_x + F \cdot \nabla_v \right]f = 0 $$

($\nabla_x$ trivially acts on $v$ and $\nabla_v$ trivially acts on $F(t,x)$) and the result follows. 

To be precise, one has to interpret $|f|^p = \lim_{\epsilon\to 0} (\sqrt{\epsilon + |f|^2})^{p}$, and approximate your $f\in L^p$ with $f \in \mathcal{S}\cap L^p$ or similar. 

For the case $p = \infty$ it suffices to notice that $f$ is constant on the integral curves of the vector field $(1, v, F(t,x))$ on $\mathbb{R}\times\mathbb{R}^n \times \mathbb{R}^n$. 

Finally, something stronger is true: let $G:\mathbb{R}\to\mathbb{R}$ be smooth, then $\int G\circ f ~dz$ is invariant in time if $f$ solves Vlasov.