In this work, we investigate the merger rate of primordial black hole–neutron star (PBH-NS) binaries in two widely studied modified gravity (MG) models: Hu–Sawicki f(R) gravity and the normal branch of Dvali–Gabadadze–Porrati gravity. In our analysis, we take into account the effects of MG on the halo properties, including halo mass function, halo concentration parameter, halo density profile, and velocity dispersion of dark matter particles. We find that these MG models, due to their stronger gravitational field induced by an effective fifth force, predict enhanced merger rates compared to general relativity. This enhancement is found to be redshift-dependent and sensitive to model parameters and PBH mass and fraction. Assuming a PBH mass range of 5–50 M
⊙, we compare the predicted merger rate of PBH-NS binaries with those inferred from LIGO–Virgo–KAGRA observations of gravitational waves (GWs). We find that the merger rates obtained from MG models will be consistent with the GW observations if the abundance of PBHs is relatively large, with the exact amount depending on the MG model and its parameter values, as well as PBH mass. We also establish upper limits on the abundance of PBHs in these MG frameworks while comparing them with the existing non-GW constraints, which can potentially impose even more stringent constraints.