As an extension of isotropic Gaussian random fields and Q-Wiener processes on d-dimensional spheres, isotropic Q-fractional Brownian motion is introduced and sample Hölder regularity in space-time is shown depending on the regularity of the spatial covariance operator Q and the Hurst parameter H. The processes are approximated by a spectral method in space for which strong and almost sure convergence are shown. The underlying sample paths of fractional Brownian motion are simulated by circulant embedding or conditionalized random midpoint displacement. Temporal accuracy and computational complexity are numerically tested, the latter matching the complexity of simulating a Q-Wiener process if allowing for a temporal error.