This paper presents a novel formula for the transition density of the Brownian motion on a sphere of any dimension and discusses an algorithm for the simulation of the increments of the spherical Brownian motion based on this formula. The formula for the density is derived from an observation that a suitably transformed radial process (with respect to the geodesic distance) can be identified as a Wright-Fisher diffusion process. Such processes satisfy a duality (a kind of symmetry) with a certain coalescent processes and this in turn yields a spectral representation of the transition density, which can be used for exact simulation of their increments using the results of Jenkins and Spanò (2017). The symmetry then yields the algorithm for the simulation of the increments of the Brownian motion on a sphere. We analyse the algorithm numerically and show that it remains stable when the time-step parameter is not too small.