The Mu<spacing diaeresis>ntz Legendre polynomials are a family of generalized orthogonal polynomials, defined by contour integral associated with a complex sequence Lambda = {lambda 0 , lambda 1, 1 , lambda 2, 2 , <middle dot> <middle dot> <middle dot> }. In this paper, we are interested in two subclasses of the Mu<spacing diaeresis>ntz Legendre polynomials. Precisely, we theoretically and numerically investigate the basic approximation properties of the Mu<spacing diaeresis>ntz Legendre polynomials for two sets of Lambda sequences: lambda k k = lambda, and lambda k k = k lambda + q for some lambda and q. 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 Mu<spacing diaeresis>ntz polynomials, a Galerkin spectral method based on the Mu<spacing diaeresis>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.