The Müntz Legendre polynomials are a family of generalized orthogonal polynomials, defined by contour integral associated with a complex sequence Λ={λ0,λ1,λ2,⋯}\Lambda =\{\lambda _{0},\lambda _{1},\lambda _{2},\cdots \}. In this paper, we are interested in two subclasses of the Müntz Legendre polynomials. Precisely, we theoretically and numerically investigate the basic approximation properties of the Müntz Legendre polynomials for two sets of Λ\Lambda sequences: λk=λ\lambda _{k}=\lambda, and λk=kλ+q\lambda _k=k\lambda +q for some λ\lambda and qq. First, the projection and interpolation errors are analyzed and numerically tested for each of the two subclasses of polynomials, and some error estimates are derived for functions in non-uniformly weighted Sobolev spaces. Then, in order to demonstrate the applicability of the Müntz polynomials, a Galerkin spectral method based on the Müntz Legendre polynomials is proposed to solve the time-space fractional differential equation. The obtained numerical results show that the proposed method leads to an exponential convergence rate even if the exact solutions are not smooth. This is opposed to low order algebraic convergence if traditional orthogonal polynomials are used.