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.