  We use separate universe simulations with the IllustrisTNG galaxy formation
model to predict the local PNG bias parameters $b_\phi$ and $b_{\phi\delta}$ of
atomic neutral hydrogen, ${\rm H_I}$. These parameters and their relation to
the linear density bias parameter $b_1$ play a key role in observational
constraints of the local PNG parameter $f_{\rm NL}$ using the ${\rm H_I}$ power
spectrum and bispectrum. Our results show that the popular calculation based on
the universality of the halo mass function overpredicts the $b_\phi(b_1)$ and
$b_{\phi\delta}(b_1)$ relations measured in the simulations. In particular, our
results show that at $z \lesssim 1$ the ${\rm H_I}$ power spectrum is more
sensitive to $f_{\rm NL}$ compared to previously thought ($b_\phi$ is more
negative), but is less sensitive at other epochs ($b_\phi$ is less positive).
We discuss how this can be explained by the competition of physical effects
such as that large-scale gravitational potentials with local PNG (i) accelerate
the conversion of hydrogen to heavy elements by star formation, (ii) enhance
the effects of baryonic feedback that eject the gas to regions more exposed to
ionizing radiation, and (iii) promote the formation of denser structures that
shield the ${\rm H_I}$ more efficiently. Our numerical results can be used to
revise existing forecast studies on $f_{\rm NL}$ using 21cm line-intensity
mapping data. Despite this first step towards predictions for the local PNG
bias parameters of ${\rm H_I}$, we emphasize that more work is needed to assess
their sensitivity on the assumed galaxy formation physics and ${\rm H_I}$
modeling strategy.
