Series expansions of isotropic Gaussian random fields on $mathbb{S}^2$ with independent Gaussian coefficients and localized basis functions are constructed. Such representations provide an alternative to the standard Karhunen-Lo`eve expansions of isotropic random fields in terms of spherical harmonics. Their multilevel localized structure of basis functions is especially useful in adaptive algorithms. The basis functions are obtained by applying the square root of the covariance operator to spherical needlets. Localization of the resulting covariance-dependent multilevel basis is shown under decay conditions on the angular power spectrum of the random field. In addition, numerical illustrations are given and an application to random elliptic PDEs on the sphere is analyzed.