We demonstrate that baryonification algorithms, which displace particles in gravity-only simulations according to physically motivated prescriptions, can simultaneously capture the impact of baryonic physics on the two and three-point statistics of matter. Specifically, we show that our implementation of a baryonification algorithm jointly fits the changes induced by baryons on the power spectrum and equilateral bispectrum on scales up to $k = 5\rm h\, {\rm Mpc}^{-1}$ and redshifts 0 ≤ z ≤ 2, as measured in six different cosmological hydrodynamical simulations. The accuracy of our fits is typically $\sim 1{{\ \rm per\ cent}}$ for the power spectrum, and for the equilateral and squeezed bispectra, which somewhat degrades to $\sim 3{{\ \rm per\ cent}}$ for simulations with extreme feedback prescriptions. Our results support the physical assumptions underlying baryonification approaches and encourage their use in interpreting weak gravitational lensing and other cosmological observables.