Estimating the locations and spatial extents of brain sources poses a long-standing challenge for electroencephalography and magnetoencephalography (E/MEG) source imaging. In the present work, a novel source imaging method, Bayesian Electromagnetic Spatio-Temporal Imaging of Extended Sources (BESTIES), which is built upon a Bayesian framework that determines the spatio-temporal smoothness of source activities in a fully data-driven fashion, is proposed to address this challenge. In particular, a Markov Random Field (MRF), which can precisely capture local cortical interactions, is employed to characterize the spatial smoothness of source activities, the temporal dynamics of which are modeled by a set of temporal basis functions (TBFs). Crucially, all of the unknowns in the MRF and TBF models are learned from the data. To accomplish model inference efficiently on high-resolution source spaces, a scalable algorithm is developed to approximate the posterior distribution of the source activities, which is based on the variational Bayesian inference and convex analysis. The performance of BESTIES is assessed using both simulated and actual human E/MEG data. Compared with L2-norm constrained methods, BESTIES is superior in reconstructing extended sources with less spatial diffusion and less localization error. By virtue of the MRF, BESTIES also overcomes the drawback of over-focal estimates in sparse constrained methods.