**Edit:** unfortunately I did botch the crucial computation. So, surprisingly, the map $f$ defined below need not be locally Lipschitz in the $C^{k,\alpha}$ space.

Let me give an example in $1$-dimension: let $v(x)=\frac23 |x|^{3/2}\in C^{1,\frac12}$, and consider for small $\varepsilon$ the functions $g(x)=\varepsilon +x$ and $h(x)=x$. Then obviously $\lVert g-h\rVert_{C^{1,\frac12}}=\varepsilon$, but we have for $x>0$
$$(v\circ g-v\circ h)'(x)=\sqrt{\varepsilon+x}-\sqrt{x}=\frac{\varepsilon}{\sqrt{\varepsilon+x}+\sqrt{x}} \to \sqrt{\varepsilon} \quad\mbox{when } x\to 0$$
so $\lVert v\circ g-v\circ h\rVert_{C^{1,\frac12}} \ge \sqrt{\varepsilon}$.

Even if this counter-example does not rule out the possibility that what I asked holds (the flow of $v$ above has the wanted property), it still makes me doubt.

**Previous version (hopefully still useful in other regularities)**

Since references seem elusive, let me propose a simple proof. The idea is simply to use the Picard–Lindelöf theorem to the right object.

First, since $M$ is compact it has positive injectivity radius, and its (global) exponential map $\exp:TM\to M$ induces a diffeomorphism $(x,u)\mapsto (x,\exp_x(u))$ from a neighborhood of the zero section of $TM$, to a neighborhood of the diagonal in $M\times M$ (of course this is already the Picard–Lindelöf theorem, but in the more classical smooth regularity and without the need for the "global derivative" part of the question). Pulling back by $\exp$, we can identify the diffeomorphisms of $M$ which are uniformly close to the identity with vector fields (which are then uniformly close to zero). We use the letter $\Phi$ to denote a diffeomorphism seen as a point the consequent open set $\Omega$ of $\Gamma^{k,\alpha}$ (the space of $C^{k,\alpha}$ vector fields), and the letter $V$ to denote an element of $\Gamma^{k,\alpha}$, seen as a tangent vector to $\Omega$.

Let $v\in \Gamma^{k,\alpha}$ be our given vector field, and define $f:\Omega\to \Gamma^{k,\alpha}$ by
$$f(\Phi) = v\circ \Phi$$
Then $f$ is locally Lipschitz in the $C^{k,\alpha}$ norm, which makes $\Gamma^{k,\alpha}$ a Banach space. This follows from a certainly classical, but somewhat tedious computation I hope I got right (the same needed to show that the $C^{k,\alpha}$ regularity is stable by product and *by composition*; this last item requires $k\ge1$)

Applying the Picard-Lindelöf theorem to the differential equation $\Phi'(t)=f(\Phi(t))$ in $\Omega\subset \Gamma^{k,\alpha}$ thus yields a unique maximal solution starting at $\Phi(0)=\mathrm{Id}_M$. This solution is obviously the flow of $v$, and now the wanted derivative in $C^{k,\alpha}$ norm follows from the equation itself.