The tetragonal to monoclinic phase transformation and its toughening effect on the Mode I crack propagation in single crystalline zirconia is investigated by a coupled phase field model. The developed phase field model can simulate the evolution of the tetragonal to monoclinic phase transformation and crack propagation simultaneously. Two tetragonal lattice orientations: 0° and 90°, are studied and compared with the case that the material is untransformable. The numerical results show that the phase transformation always initiates from the crack tip. The phase transformation dramatically reduces the normal stress near the crack tip and slows down the crack growth speed for both orientation angles. A closure of the crack tip may even be observed shortly after the phase transformation, which is caused by the transformation stain and the stress re-distribution at the junction of the monoclinic phase and the crack phase. The lattice orientations have great influences on the final monoclinic microstructure and the transformation toughening effect. This work shows the validity of the phase field method to simulate the microstructure evolution during the crack propagation and it can serve as a robust tool to investigate the transformation toughening behavior of Zirconia ceramics.