Linear magnon-phonon coupling hybridizes magnon and phonon bands at the same energy and momentum, resulting in an anticrossing signature. This hybrid quasiparticle benefits from a long phonon lifetime and efficient magnon transport, showing great potential for spintronics and quantum information science applications. In this paper, we present an efficient and accurate first-principles approach for calculating linear magnon-phonon couplings. We first calculate the magnon spectra from linear spin-wave theory with spin Hamiltonian and firstprinciples exchange constants, which compared well with time-dependent density-functional theory. We then obtain the magnon-phonon coupling from the derivative of off-diagonal exchange constants in real space, calculated from the Hellmann-Feynman forces of the spin-constrained configurations, avoiding the use of cumbersome finite-difference methods. Our implementation allows calculating coupling coefficients at an arbitrary wave vector in the Brillouin zone in a single step, through Fourier interpolation of real-space supercell calculations. We verify our implementation through two-dimensional magnetic systems, monolayer CrI3, in agreement with experiments, and extend its application to monolayer CrTe2. We emphasize the role of nonmagnetic atoms in superexchange interactions and magnon-phonon coupling, which have been overlooked previously. We suggest effective tuning of magnon-phonon coupling through strain, doping, and terahertz excitations, for spintronics and quantum magnonics applications.