The nuclear-electronic orbital (NEO) method is combined with vibronic coupling theory to calculate hydrogen tunneling splittings in polyatomic molecules. In this NEO-vibronic coupling approach, the transferring proton and all electrons are treated quantum mechanically at the NEO level, and the other nuclei are treated quantum mechanically using vibronic coupling theory. The dynamics of the molecule are described by a vibronic Hamiltonian in a diabatic basis of two localized nuclear-electronic states for the electrons and transferring proton. This ab initio approach is computationally practical and efficient for relatively large molecules, and the accuracy can be improved systematically. The NEO-vibronic coupling approach is used to calculate the hydrogen tunneling splitting for malonaldehyde. The calculated tunneling splitting of 24.5 cm(-1) is in excellent agreement with the experimental value of 21.6 cm(-1). This approach also enables the identification of the dominant modes coupled to the transferring hydrogen motion and provides insight into their roles in the hydrogen tunneling process.