Is every geodesic-preserving diffeomorphism an isometry?

Let $M$ be a closed $n$-dimensional Riemannian manifold.

Let $f:M \to M$ be a diffeomorphism and suppose that for every (parametrized) geodesic $\gamma$, $f \circ \gamma$ is also a (parametrized) geodesic.

Must $f$ be an isometry?

An equivalent condition on $f$ is that $\nabla df=0$ where $\nabla=\nabla^{T^*M} $ $ \otimes \nabla^{f^*TM}$ is the relevant tensor product connection.

Note that this equivalent assumption implies that $df$ has constant singular values, and in particular that the Jacobian $\det(df)$ is constant, hence it must be $1$ (since $f$ was assumed to be a diffeomorphism.) Thus $f$ is volume-preserving.

This can certainly be false for manifolds with nonempty boundary in general, as the following example shows:

Let $0<a<b$, and set $ M=D_{a,b}=\biggl\{(x,y) \,\biggm | \, \frac{x^2}{a^2} + \frac{y^2}{b^2} \le 1 \biggr\} $ to be the ellipse with diameters $a,b$, endowed with the standard Euclidean metric (induced by $\mathbb{R}^2$). Then there exists $A \in \operatorname{SL}_2(\mathbb R) \setminus \operatorname{SO}(2)$ such that $AD_{a,b}=D_{a,b}$, and $A$ clearly preserves geodesics (it maps straight lines to straight lines.)

Indeed one can take $A$ to be of the form $$ A =A_{\theta}:= \begin{pmatrix} a& 0 \\ 0 & b \end{pmatrix} \begin{pmatrix} \cos\theta & -\sin\theta \\ \sin\theta & \cos \theta \end{pmatrix}\begin{pmatrix} 1/a& 0 \\ 0 & 1/b \end{pmatrix}= \begin{pmatrix} \cos\theta & -\frac ab \sin\theta \\ \frac ba \sin\theta & \cos \theta \end{pmatrix}. $$

For $M=\mathbb{S}^n $ the answer is positive, by this answer here.

No. As a counterexample, take the torus $\mathbb{T}^n=\mathbb{R}^n/\mathbb{Z}^n$ where $\mathbb{R}^n$ is endowed with the Euclidean metric. Since linear maps preserve geodesics in $\mathbb{R}^n$, $GL(n,\mathbb{Z})$ acts by geodesic-preserving diffeomorphisms on $\mathbb{T}^n$, and many of these are not isometries. For example, consider the map $\mathbb{T}^2\to\mathbb{T}^2$ defined by $(x,y)\mapsto(x,x+y)$ modulo $1$.

Qualitatively speaking, the requirements on $f$ do not constrain the shear components of $df$, though I suppose there may be global or even local obstructions to the existence of such shear maps.