We use hydrodynamical separate universe simulations with the IllustrisTNG model to predict the local primordial non-Gaussianity (PNG) bias parameters b φ and b φδ , which enter at leading order in the galaxy power spectrum and bispectrum. This is the first time that b φδ is measured from either gravity-only or galaxy formation simulations. For dark matter halos, the popular assumption of universality overpredicts the b φδ (b 1 ) relation in the range 1 b 1 3 by up to ∆b φδ ∼ 3 (b 1 is the linear density bias). The adequacy of the universality relation is worse for the simulated galaxies, with the relations b φ (b 1 ) and b φδ (b 1 ) being generically redshift-dependent and very sensitive to how galaxies are selected (we test total, stellar and black hole mass, black hole mass accretion rate and color). The uncertainties on b φ and b φδ have a direct, often overlooked impact on the constraints of the local PNG parameter f nl , which we study and discuss. For a survey with V = 100Gpc 3 /h 3 at z = 1, uncertainties ∆b φ 1 and ∆b φδ 5 around values close to the fiducial can yield relatively unbiased constraints on f nl using power spectrum and bispectrum data. We also show why priors on galaxy bias are useful even in analyses that fit for products f nl b φ and f nl b φδ . The strategies we discuss to deal with galaxy bias uncertainties can be straightforwardly implemented in existing f nl constraint analyses (we provide fits for some of the bias relations). Our results motivate more works with galaxy formation simulations to refine our understanding of b φ and b φδ towards improved constraints on f nl .