Elliptical subsurface cracks are one of the probable types of cracks that occur in engineering structures, especially under rolling contact fatigue. Due to the non-symmetrical geometry, coupling of the fracture modes occurs in an elliptical subsurface crack and the crack under shear loading will experience all fracture modes. This paper investigates the coupling of the fracture modes of elliptical subsurface cracks under uniform shear loadings in two directions. First, a three-dimensional parametric finite element model of a crack in an infinite space has been developed and validated. Then, by moving the crack close to surface, mixed mode stress intensity factors (SIFs) have been calculated for cracks with aspect ratios of 0.2-1.0 and ratios of crack depth to crack length of 0.05-1.0. Based on the results, coupling of the fracture modes occur considerably when the crack depth becomes less than crack length. By decreasing of the crack depth from infinite to 0.05, shear SIFs and KImax/KIImax ratio increase at least up to 65% and 90%, respectively. Six equations for SIFs of the subsurface cracks under uniform shear loadings in two directions have been obtained by fitting to the finite element results. These equations can be used efficiently in high accurate calculation of the SIFs for subsurface cracks with any aspect ratio and depth under uniform shear loading with any direction.