Let $(M,h)$ be a compact connected smooth Riemannian manifold with $\partial M=\emptyset$ and $h$ smooth metric. Let $p_t(x,y)$ designate its heat kernel $(t>0,x,y\in M)$.
In the Euclidean case, one would have that $$ \nabla_x p_t(x,y) = -\nabla_y p_t(x,y)\,. $$
My question is: does a similar identity hold on $M$, through parallel transport (since the left hand side belongs to $T_xM$ and the right hand side to $T_yM$)?
(I am aware of Bismut formula which enables one to represent $\nabla_x p_t(x,y)$ in a probabilistic setting, but I don't know how to proceed from it.)
Certainly for symmetric spaces we have the reflection isometry $x\leftrightarrow y$ so it holds just like $\mathbb{R}^n$.
But I don't believe it works in general. Demanding the symmetry $\nabla_x p_t(x,y)=-P_{y\to x}\nabla_y p_t(x,y)$ along any geodesic joining $x$ and $y$ is basically equivalent to $$ P_{\gamma(s)\to\gamma(0)}\frac{D}{ds}\log\phi(\gamma(s)) $$ being an odd function in $s$ for every geodesic $\gamma$, and every eigenfunction $\phi$ of the Laplacian, which looks like too much to ask for.