Functional maps are efficient representations of shape correspondences, that provide matching of real-valued functions between pairs of shapes. Functional maps can be modelled as elements of the Lie group $SO(n)$ for nearly isometric shapes. Synchronization can subsequently be employed to enforce cycle consistency between functional maps computed on a set of shapes, hereby enhancing the accuracy of the individual maps. There is an interest in developing synchronization methods that respect the geometric structure of $SO(n)$, while introducing a probabilistic framework to quantify the uncertainty associated with the synchronization results. This paper introduces a Bayesian probabilistic inference framework on $SO(n)$ for Riemannian synchronization of functional maps, performs a maximum-a-posteriori estimation of functional maps through synchronization and further deploys a Riemannian Markov-Chain Monte Carlo sampler for uncertainty quantification. Our experiments demonstrate that constraining the synchronization on the Riemannian manifold $SO(n)$ improves the estimation of the functional maps, while our Riemannian MCMC sampler provides for the first time an uncertainty quantification of the results.