Let $M$ be a compact Riemannian manifold and let $\mathrm{d}\mathbb{W}^{yx;T}(\gamma)$ denote the Brownian Bridge measure, i.e. the Wiener measure on the paths that travel from $x$ to $y$ in time $T$ (in such a way that the integral over the function $1$ gives back the heat kernel $e^{-T\Delta}(x, y)$ on $M$).
Now if $f$ and $g$ are some smooth functions on $M$, then we have $$ \int f\bigl(\gamma(s_2)\bigr)\, g\bigl(\gamma(s_1)\bigr)\, \mathrm{d}\mathbb{W}^{yx;T}(\gamma) = \begin{cases} \bigl(e^{-(T-s_2)\Delta} M_f e^{-(s_2-s_1)\Delta}M_g e^{-s_1\Delta}\bigr)(x, y) & \text{if}~~0 \leq s_1 \leq s_2 \leq T\\ \bigl(e^{-(T-s_1)\Delta} M_g e^{-(s_1-s_2)\Delta}M_f e^{-s_2\Delta}\bigr)(x, y) & \text{if}~~ 0 \leq s_2 \leq s_1 \leq T,\end{cases}$$ where $M_f$ and $M_g$ denote the multiplication operators that multiply with $f$, respectively with $g$. That is, the path integral has a built-in "time-ordering" operator, in physicist's slang.
Now if we additionally have some metric vector bundle over $M$ with a metric connection $\nabla$, we can consider the operator $L = \nabla^* \nabla$, whose heat kernel is now given by the Feynman-Kac-Ito formula $$ e^{-TL}(x, y) = \int [\gamma\|_0^T]\, \mathrm{d}\mathbb{W}^{yx;T}(\gamma),$$ where $[\gamma\|_0^T]$ denotes stochastic parallel transport along the path $\gamma$.
Similar to the above, if $F$ and $G$ are given endomorphism fields on $M$, we can consider the integral $$\int [\gamma\|_{s_2}^T]\,F\bigl(\gamma(s_2)\bigr)\,[\gamma\|_{s_1}^{s_2}]\,G\bigl(\gamma(s_1)\bigr)\,[\gamma\|_0^{s_1}]\, \mathrm{d}\mathbb{W}^{yx;T}(\gamma).$$ If $0 \leq s_1 \leq s_2 \leq T$, the integral is easily found to equal the kernel of the operator $e^{-(T-s_2)L}M_F e^{-(s_2-s_1)L}M_Ge^{-s_1L}$, evaluated at $(x, y)$. However, I couldn't see what happens in the case that $0 \leq s_2 \leq s_1 \leq T$. For example, it is not so easy to pull the integrals apart, because the integrands are not independent.
