For notational convenience I'm going to write $v(x)$ instead of $\vec{xy},$ which defines a vector field $v \in \Gamma(TM)$. I will use $D$ to denote the Levi-Civita connection and reserve $\nabla$ for the gradient. When I write $A \lesssim B$ it means $A \le C\, B$ for some constant $C$ depending only on the dimension.
Let $\sigma : [0,1] \to M$ be the geodesic segment joining $x$ to $x'$ and $V\in\Gamma(\sigma^* TM)$ the vector field along $\sigma$ obtained by parallel transporting $v(x) \in T_{x}M$. We are trying to estimate $E=|V(1) - v(\sigma(1))|.$
We can do this by studying the cumulative error $f(t) = |V(t) - v(t)|,$ which clearly satisfies $f(0) = 0$ and $f(1) = E.$ (I'm now using $v(t)$ to denote the restriction $v \circ \sigma(t).$) Remembering that $V$ is parallel, we can differentiate this using metric-compatibility to get $$\frac{df}{dt} = \frac 1{f}\left\langle V - v,- \frac{Dv}{dt}\right\rangle,$$ which we can estimate using Cauchy-Schwarz to get $f'(t) \le |Dv/dt|$ and thus $$E = f(0) + \int_0^1 f'\le\sup \left|\frac{Dv}{dt}\right| \le d(x,x') \sup_{\sigma} |Dv|,$$ where the last step comes from the fact that $|\dot\sigma|=d(x,x').$ Thus we just need to estimate the covariant derivative $Dv.$
Letting $r(p) = d(y,p)$ denote the distance to $y$, note that $v = - r \nabla r$ since $r$ has unit-length gradient pointing along the geodesic rays from $y.$ Estimating the derivative of $v$ thus requires an estimate on the Hessian of the distance function $H = D^2r$. More precisely, differentiating $v = -r \nabla r$ yields $Dv = -dr \otimes \nabla r-r H^\sharp,$ so for $p$ in the image of $\sigma$ we have $$|Dv(p)| \lesssim 1 + \max(d(x,y),d(x',y)) |H(p)|.$$ Assuming the sectional curvature of $M$ is bounded below by some constant $K$, the Hessian Comparison Theorem (see e.g. Hamilton's Ricci Flow by Chow et al) tells us that $|H(p)| \lesssim |H_K(r(p))|$ where $H_K(r)$ is the mean curvature of the sphere of radius $r$ in the space of constant curvature $K$. For example, in the worst case $K < 0$ we have $$H_K(r) = (n-1)\sqrt{|K|}\coth\left(\sqrt {|K|} r\right).$$ Since $H_K$ is monotone, we thus have an estimate $$E\le C(n)\Big(1+\max(d,d')\max(|H_K(d)|,|H_K(d')|)\Big)d(x,x')$$ where $d = d(x,y), d' = d(x',y).$