Let $M,N$ be smooth $d$-dimensional Riemannian manifolds.
Suppose $f:M \to N$ is a $C^1$ isometry (i.e $f$ is continuously differentiable, and $df_p$ is an isometry for every $p \in M$).
I think that $f$ must be smooth (I describe a way of showing this below). Is there a simpler way?
It suffices to show $f$ maps geodesic to geodesic, that is: $\alpha$ is a geodesic $\Rightarrow$ $f \circ \alpha$ is a geodesic.
Indeed, if this is the case and $\alpha(0)=p,\alpha'(0)=v$, then $f \circ \alpha(0)=f(p),(f \circ \alpha)'(0)=df_p(v)$, so by uniqueness of geodeics, we conclude that $f \circ exp_p^M(tv)=f \circ \alpha(t)=exp_{f(p)}^N(t \cdot df_p(v))$. In other words, $$f \circ exp_p^M=exp_{f(p)}^N\circ df_p,$$
so $f$ is a local diffeomorphism, hence smooth.
The problem is: How do we know $f$ maps geodesic to geodesic?
If we take the differntial definition, $\nabla \dot\alpha=0$, then we hit a problem since a-priori $f \circ \alpha$ is differentiable only once (the geodesic equation is of second order).
However, we could use the characterization of geodesics as length minimizers, but this requires that we know that geodesics "beat" any $C^1$ paths, not just (piecewise) smooth ones. This is essentially saying that the Riemannian distance does not change when using only $C^1$ paths.
Another way to proceed is to use the Myers-Steenrod theorem:
We consider the lenght structure $(\mathcal{A}_M,L_M)$ where $\mathcal{A}$ is the set of all $C^1$ paths in $M$, and $L_M:\mathcal{A}_M \to \mathbb{R}$ is the standard Riemannian length: $L(\gamma)=\int \|\dot \gamma(t)\|$. (and similarly for $N$).
By our assumption, $f$ maps paths in $\mathcal{A}_M$ to paths in $\mathcal{A}_N$, and it's also length-preserving, i.e an arcwise isometry.
By the inverse function theorem, $f$ is a local homeomorphism, hence a local isometry w.r.t the intrinsic distances defined using the class of $C^1$ paths. (By intrinsic I mean we use only paths that stay inside the neighbourhood and take the associated distance, not the restricted distance from the manifolds).
Since the Riemannian distance does not change when using only $C^1$ paths, it follows that $f$ is a local isometry w.r.t the intrinsic distances (defined as usual using the class of piecsewise smooth paths paths).
Now the smoothness of $f$ follows from the Myers-Steenrod theorem.
As mentioned in the quesiton, it's enough to show that if $\alpha$ is a geodesic, then $f \circ \alpha$ is a geodesic.
Let $p \in M$, and let $\alpha$ be a geodesic emanating from $p$.
Then for sufficiently small $t$, it holds that $L(\alpha|_{[0,t]}) \le L(\beta|_{[0,t]})$ for any $C^1$ path $\beta$ connecting $\alpha(0),\alpha(t)$. (This follows from the standard proof that geodesics are locally length minimizing, which works for $C^1$ paths).
By the inverse funciton theorem there exist open neighbourhoods $p \in U,f(p) \in V$ s.t $\, \, f:U \to V$ is invertible.
We claim that $f \circ \alpha$ is also length minimizing for sufficiently small $t$. Indeed, if there was a shorter $C^1$ path $c$ between $f(p),f(\alpha(t))$, we could assume W.L.O.G that it lies completely inside $V$, hence $f^{-1} \circ c$ was a $C^1$ path shorter than $\alpha$, a contradiction.
Thus, $f \circ \alpha$ is a geodesic, as required. (Since it's length minimizing and parametrized by arc-length).