  We use hydrodynamical separate universe simulations with the IllustrisTNG
model to predict the local primordial non-Gaussianity (PNG) bias parameters
$b_{\phi}$ and $b_{\phi\delta}$, which enter at leading order in the galaxy
power spectrum and bispectrum. This is the first time that $b_{\phi\delta}$ is
measured from either gravity-only or galaxy formation simulations. For dark
matter halos, the popular assumption of universality overpredicts the
$b_{\phi\delta}(b_1)$ relation in the range $1 \lesssim b_1 \lesssim 3$ by up
to $\Delta b_{\phi\delta} \sim 3$ ($b_1$ is the linear density bias). The
adequacy of the universality relation is worse for the simulated galaxies, with
the relations $b_{\phi}(b_1)$ and $b_{\phi\delta}(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_{\phi}$ and $b_{\phi\delta}$ have a direct, often
overlooked impact on the constraints of the local PNG parameter $f_{\rm NL}$,
which we study and discuss. For a survey with $V = 100{\rm Gpc}^3/h^3$ at
$z=1$, uncertainties $\Delta b_{\phi} \lesssim 1$ and $\Delta b_{\phi\delta}
\lesssim 5$ around values close to the fiducial can yield relatively unbiased
constraints on $f_{\rm 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_{\rm NL} b_{\phi}$ and $f_{\rm NL} b_{\phi\delta}$. The strategies
we discuss to deal with galaxy bias uncertainties can be straightforwardly
implemented in existing $f_{\rm 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_{\phi}$ and
$b_{\phi\delta}$ towards improved constraints on $f_{\rm NL}$.
