To obtain the result $\frac{\text d \rho }{\text d t}=0$ you need two facts: the first is that the hamiltonian flow preserves the volume of phase space. The second fact is the conservation of probability, that is, the probability that the system is found in a volume $U$ at time $t=0$ equals the probability of finding it within $\Phi _t U$ at time $t$, where $\Phi _t$ denotes the hamiltonian flow. This is a direct consequence of the deterministic nature of classical mechanics: the two propositions “$(p(0),q(0))\in U$” and “$(p(t),q(t))\in \Phi _t U$” are equivalent.
Using conservation of probability, for an arbitrary volume $U$ we can write an equation: $$\int _U \rho(p,q,0) \text d p \text d q=\int _{\Phi _t U} \rho(p,q,t)\text dp \text d q .$$
By Jacobi's theorem: $$\int _{\Phi_t U} \rho (p,q,t)\text d p \text d q=\int _U\rho (\Phi _t (p,q),t)\text J_{\Phi _t}d p \text d q.$$
The Jacobian $J_{\Phi _t}=1$, because the flow preserves volumes. It follows that: $$\int _U \rho (p,q,0)\text d p \text d q =\int _U \rho (\Phi _t (p,q),t)\text d p \text d q,$$
and, since the volume $U$ was arbitrary, $\rho (p,q,0)=\rho (\Phi _t (p,q),t)$, or $\text d\rho /\text d t=0$.
To answer the response in the comment, yes this is merely a coincidence. The Hamiltonian mechanics theorem is a "dynamical statement". It says if you have a Hamiltonian $H$ and look at the integral flow $\Phi_{H,t}$, then this flow leaves the symplectic form invariant $\Phi_{H,t}^*\omega=\omega$ for all $t$. It is dynamical in the sense that it talks about the Hamiltonian $H$ itself, its vector field $X_H:=\omega^{\sharp}(dH)$ and the corresponding integral flow. Meanwhile, the symplectic form $\omega$ talks about the geometry of phase space. So, this version of Liouville's theorem says roughly that the geometry of phase space is unaffected by the dynamics.
The one about complex analysis is just an emphasis of the rigid behavior of holomorphic (or even harmonic) functions. There are so many theorems in math/physics which deal with some quantity being constant/ being conserved under certain actions. Some have relationships with one another, but not always (and this isn't one of those cases). The fact that bounded entire functions are constant is a special case of the following result (whose proof is very similar to the usual proof of Liouville' theorem).
Let $f:\Bbb{C}\to\Bbb{C}$ be an entire function, and suppose it satisfies an estimate of the form $|f(z)|\leq C|z|^{\alpha}$, for some $C,\alpha\geq 0$. Then, $f$ must be a polynomial of degree at most $\lfloor{\alpha}\rfloor$.
This says if $f$ satisfies some kind of bound, then $f$ has to be a polynomial of not-too large degree. In particular if you take $0\leq \alpha<1$, you can deduce $f$ is constant (the special case $\alpha=0$ corresponds to the usual Liouville theorem). As you can see this theorem is more of a restriction on how crazy holomorphic functions can be (namely, not very crazy at all... compared to smooth functions on $\Bbb{R}$ which can have crazy behavior).
Best Answer
Here's something I believe is a simple proof. Unfortunately it uses a little bit of cohomology.
Consider the canonical 2-form in extended phase space $T^*M \times \mathbb{R}$
$$\omega = \sum_{i=1}^N dq_i \wedge dp_i - dH(\vec{q},\vec{p},t) \wedge dt ,$$
where $N = dim(M)$. A function $f: M \to M$ is said to be a canonical transformation iff $f^* \omega = \omega$. Thus, for any canonical transformation,
$$\sum_{i=1}^N dq_i \wedge dp_i - dH(\vec{q},\vec{p},t) \wedge dt = \sum_{i=1}^N dQ_i \wedge dP_i - dK(\vec{Q},\vec{P},T) \wedge dT ,$$
where we defined $(\vec{Q},\vec{P},T) = f^*(\vec{q},\vec{p},t)$. This means that
$$\sum_{i=1}^N q_i dp_i - H(\vec{q},\vec{p},t) dt - \left( \sum_{i=1}^N Q_i dP_i - dK(\vec{Q},\vec{P},T) dT \right) = dG ,$$
i.e, that the tautological forms associated to the canonical 2-forms are elements of the same de Rham cohomology class. By the homotopic invariance of the de Rham cohomology,
$$\oint_\gamma \left( \sum_{i=1}^N q_i dp_i - H(\vec{q},\vec{p},t) dt \right) = \oint_{\Gamma} \left( \sum_{i=1}^N Q_i dP_i - dK(\vec{Q},\vec{P},T) dT \right) ,$$
where $\Gamma$ is the image of the curve $\gamma$ by $f$.
As movement can be considered to be a canonical transformation, we just proved that the hamiltonian evolution in extended phase space needs to relate only homotopic curves. A corollary is that if you start with a compact simply-connected set in phase space and track its evolution in time, the curve that bounds this set can only evolve to others homotopic to it. As a simply connected set boundary is never homotopic to a non-simply connected one, we just showed that hamiltonian evolution takes simply-connected sets to simply-coonected sets.
I'm really sorry if this is too technical. I'm too ignorant to offer a proof with simpler arguments (which points out how limited is my knowledge on the subject).
EDIT: Small correction on dimension.