Isometries preserve geodesics Let $f$ be an isometry (i.e a diffeomorphism which preserves the Riemannian metrics) between Riemannian manifolds $(M,g)$ and $(N,h).$ 
One can argue that $f$ also preserves the induced metrics $d_1, d_2$ on $M, N$ from $g, h$ resp. that is, $d_1(x,y)=d_2(f(x),f(y))$ for $x,y \in M.$ Then, it's easy to show that $f$ sends geodesics on $M$ to geodesics on $N,$ using the length minimizing property of geodesics and that $f$ is distance-preserving. 
My question, 

Is it possible to derive the result without using the distance-preserving property of isometries, by merely the definition?

What I have found so far;
Let $\gamma : I \to M$ be a geodesic on $M$, i.e. $\frac{D}{dt}(\frac{d\gamma}{dt})=0,$ where $\frac{D}{dt}$ is the covariant derivative and $ \frac{d\gamma}{dt}:=d\gamma(\frac{d}{dt}).$ Let $t_0 \in I,$ we have to show that $\frac{D}{dt}(\frac{d(f \circ \gamma)}{dt})=0$ at $t=t_0,$ or $\frac{D}{dt}(df_{\gamma(t_0)}(\frac{d \gamma}{dt}|_{t=t_0}))=0.$
We also know that 
$$ \langle \frac{d \gamma}{dt}|_{t=t_0},\frac{d \gamma}{dt}|_{t=t_0}\rangle_{\gamma(t_0)}= \langle df_{\gamma(t_0)}(\frac{d \gamma}{dt}|_{t=t_0}), df_{\gamma(t_0)}(\frac{d \gamma}{dt}|_{t=t_0})\rangle_{f(\gamma(t_0))}.$$
Since $\frac{d}{dt} \langle \frac{d \gamma}{dt}|_{t=t_0},\frac{d \gamma}{dt}|_{t=t_0}\rangle_{\gamma(t_0)}=2 \langle \frac{D}{dt}(\frac{d \gamma}{dt}|_{t=t_0}),\frac{d \gamma}{dt}|_{t=t_0}\rangle_{\gamma(t_0)}=0,$ therefore
$$\frac{d}{dt} \langle df_{\gamma(t_0)}(\frac{d \gamma}{dt}|_{t=t_0}), df_{\gamma(t_0)}(\frac{d \gamma}{dt}|_{t=t_0})\rangle_{f(\gamma(t_0))}=$$
$$2\langle \frac{D}{dt}(df_{\gamma(t_0)}(\frac{d \gamma}{dt}|_{t=t_0})), df_{\gamma(t_0)}(\frac{d \gamma}{dt}|_{t=t_0})\rangle _{f(\gamma(t_0))}=0.$$
How can I conclude from $\langle\frac{D}{dt}(df_{\gamma(t_0)}(\frac{d \gamma}{dt}|_{t=t_0})), df_{\gamma(t_0)}(\frac{d \gamma}{dt}|_{t=t_0}) \rangle _{f(\gamma(t_0))}=0$ that $\frac{D}{dt}(df_{\gamma(t_0)}(\frac{d \gamma}{dt}|_{t=t_0}))=0?$
 A: It is actually an exercise in the Lee's book, I try to do it by following the hint. First, you have to understand the naturality of Riemannian connection, then everything will be clear. I like using $\nabla_{\frac{d}{dt}}$ instead of $D_t$ here.
First, Define an operator $\varphi^*\tilde{\nabla}_{\frac{d}{dt}}:\mathcal{J}(\gamma)\rightarrow\mathcal{J}(\gamma)$ by $(\varphi^*\tilde{\nabla}_{\frac{d}{dt}})V=\varphi^*(\tilde{\nabla}_{\frac{d}{dt}}(\varphi_*V))$. Suppose $V$ is a vector field along $\gamma.$ Show this operator satisfies three properties in the Lemma 4.9. in Lee's Riemannian Geometry.
(1) Linearity over $\mathcal{R}$:$$(\varphi^*\tilde{\nabla}_{\frac{d}{dt}})(aV+bW)=a(\varphi^*\tilde{\nabla}_{\frac{d}{dt}})V+b(\varphi^*\tilde{\nabla}_{\frac{d}{dt}})W$$ (2) Product rule: $$(\varphi^*\tilde{\nabla}_{\frac{d}{dt}})(fV)=\frac{df}{dt}V+f(\varphi^*\tilde{\nabla}_{\frac{d}{dt}})V.$$
  (3) If $V$ is extendible, then for any extension $\tilde{V}$ of $V$,
  $$(\varphi^*\tilde{\nabla}_{\frac{d}{dt}})V=(\varphi^*\tilde{\nabla}_{\dot{\gamma}(t)})\tilde{V}.$$
  Then by uniqueness, the operator we defined above is just the unique operator $\nabla_{\frac{d}{dt}}$. that is
  $$\nabla_{\frac{d}{dt}}V=\varphi^*(\tilde{\nabla}_{\frac{d}{dt}}(\varphi_*V)) \ \text{or} \ \ \varphi_*(\nabla_{\frac{d}{dt}}V)=\tilde{\nabla}_{\frac{d}{dt}}(\varphi_*V).$$
Now If $\gamma$ is the geodesic in $M$ with initial $p$ and initial velocity $V$, i.e., $\nabla_{\frac{d}{dt}}\dot{\gamma}(t)=0$. Obviously $\varphi\circ\gamma$ is a curve in $\tilde{M}$ with initial point $\varphi(p)$ and initial velocity $\varphi_*V$, moreover, it is also the geodesic since
 $$\tilde{\nabla}_{\frac{d}{dt}}\dot{\varphi}(\gamma(t))=\tilde{\nabla}_{\frac{d}{dt}}\varphi_*\dot{\gamma}(t)=\varphi_*(\nabla_{\frac{d}{dt}}\dot{\gamma}(t))=0.$$
A: Your calculation looks like an attempt to prove the naturality of the Levi-Civita connection, the fact that @Zhen Lin implicitly points to. In the settings of the question it can be stated as
$$
\nabla^g_X{Y}=f^* \left( \nabla^{(f^{-1})^* g}_{\operatorname{d}f(X)} \operatorname{d}f(Y) \right)
$$
Notice also that in fact you are using two different connections: one for vector fields along $\gamma \colon I \rightarrow M$ induced from $\nabla^g$ on $M$, and another one for vector fields along $f \circ \gamma \colon I \rightarrow N$ induced from $\nabla^h$ on $N$. Due to the naturality property they agree, and it may be helpful to distinguish $D^g_t:=\nabla^g_{\frac{d}{dt}{\gamma}}$ and $D^h_t:=\nabla^h_{\frac{d}{dt}(f \circ \gamma)}$ in the present calculation. Indeed, using
$$\frac{\operatorname{d}}{\operatorname{d}t}(f\circ \gamma)=\operatorname{d}f(\frac{\operatorname{d}}{\operatorname{d}t}\gamma)
$$
we get
$$
D^h_t{\frac{\operatorname{d}(f\circ\gamma)}{\operatorname{d}t}} = f_* \left( D^g_t{\frac{\operatorname{d}\gamma}{\operatorname{d}t}} \right) =0
$$
