We present the first implementation of spin−orbit coupling effects in fully internally contracted second-order quasidegenerate N-electron valence perturbation theory (SO-QD-NEVPT2). The SO-QDNEVPT2 approach enables the computations of ground-and excited-state energies and oscillator strengths combining the description of static electron correlation with an efficient treatment of dynamic correlation and spin−orbit coupling. In addition to SO-QDNEVPT2 with the full description of one-and two-body spin−orbit interactions at the level of two-component Breit−Pauli Hamiltonian, our implementation also features a simplified approach that takes advantage of spin−orbit mean-field approximation (SOMF-QDNEVPT2). The accuracy of these methods is tested for the group 14 and 16 hydrides, 3d and 4d transition metal ions, and two actinide dioxides (neptunyl and plutonyl dications). The zero-field splittings of group 14 and 16 molecules computed using SO-QDNEVPT2 and SOMF-QDNEVPT2 are in good agreement with the available experimental data. For the 3d transition metal ions, the SO-QDNEVPT2 method is significantly more accurate than SOMF-QDNEVPT2, while no substantial difference in the performance of two methods is observed for the 4d ions. Finally, we demonstrate that for the actinide dioxides the results of SO-QDNEVPT2 and SOMF-QDNEVPT2 are in good agreement with the data from previous theoretical studies of these systems. Overall, our results demonstrate that SO-QDNEVPT2 and SOMF-QDNEVPT2 are promising multireference methods for treating spin−orbit coupling with a relatively low computational cost.