Case-hardening is widely used to enhance gear loading capacity. Simulation of the material gradient properties and contact characteristics are the key issues in contact fatigue analysis of case-hardened gears. In this work, a plasto-elastohydrodynamic lubrication (PEHL) model incorporating the hardness gradient and surface roughness is developed to investigate the contact performance of case-hardened gears. The generalized Reynolds equation is solved to determine film thickness and contact pressure. The plastic deformation and residual stress are obtained via the half-space eigenstrain problem solving. The Dang Van multiaxial fatigue criterion and the Euler transformation are employed to evaluate the contact fatigue parameter based on the predetermined stress field. The discrete convolution and fast Fourier transform (DC-FFT) algorithm is used for accelerating the computation. The influences of effective case depth, surface hardness and surface roughness on the contact performance are investigated. Numerical results indicate that as the surface hardness increases, the probability of fatigue crack nucleation decreases, and the depth of the crack initiation site increases. For a lower surface roughness case, the maximum von Mises stress and equivalent plastic strain appear at a deeper layer. As the surface roughness increases, the maximum values of pressure and stress increase sharply and move closer to the surface.