This paper proposes a novel framework for manifold-valued regression and establishes its consistency as well as its contraction rate. It assumes a predictor with values in the interval $[0,1]$ and response with values in a compact Riemannian manifold $M$. This setting is useful for applications such as modeling dynamic scenes or shape deformations, where the visual scene or the deformed objects can be modeled by a manifold. The proposed framework is nonparametric and uses the heat kernel (and its associated Brownian motion) on manifolds as an averaging procedure. It directly generalizes the use of the Gaussian kernel (as a natural model of additive noise) in vector-valued regression problems. In order to avoid explicit dependence on estimates of the heat kernel, we follow a Bayesian setting, where Brownian motion on $M$ induces a prior distribution on the space of continuous functions $C([0,1], M)$. For the case of discretized Brownian motion, we establish the consistency of the posterior distribution in terms of the $L_{q}$ distances for any $1 \leq q < \infty$. Most importantly, we establish contraction rate of order $O(n^{-1/4+ε})$ for any fixed $ε>0$, where $n$ is the number of observations. For the continuous Brownian motion we establish weak consistency.