Abstract We measure power spectrum response functions in the presence of baryonic physical processes using separate universe simulations with the IllustrisTNG galaxy formation model. The response functions describe how the small-scale power spectrum reacts to long-wavelength perturbations and they can be efficiently measured with the separate universe technique by absorbing the effects of the long modes into a modified cosmology. Specifically, we focus on the total first-order matter power spectrum response to an isotropic density fluctuation R1(k, z), which is fully determined by the logarithmic derivative of the non-linear matter power spectrum dlnPm(k, z)/dlnk and the growth-only response function G1(k, z). We find that G1(k, z) is not affected by the baryonic physical processes in the simulations at redshifts z < 3 and on all scales probed (k ≲ 15 h Mpc−1; i.e. length scales $\gtrsim 0.4\, {\rm Mpc}\,h^{-1}$). In practice, this implies that the power spectrum fully specifies the baryonic dependence of its response function. Assuming an idealized lensing survey set-up, we evaluate numerically the baryonic impact on the squeezed-lensing bispectrum and the lensing supersample power spectrum covariance, which are given in terms of responses. Our results show that these higher-order lensing statistics can display varying levels of sensitivity to baryonic effects compared to the power spectrum, with the squeezed bispectrum being the least sensitive. We also show that ignoring baryonic effects on lensing covariances slightly overestimates the error budget (and is therefore conservative from the point of view of parameter error bars) and likely has negligible impact on parameter biases in inference analyses.