  Star forming galaxies emit GeV- and TeV-gamma rays that are thought to
originate from hadronic interactions of cosmic-ray (CR) nuclei with the
interstellar medium. To understand the emission, we have used the moving mesh
code Arepo to perform magneto-hydrodynamical galaxy formation simulations with
self-consistent CR physics. Our galaxy models exhibit a first burst of star
formation that injects CRs at supernovae. Once CRs have sufficiently
accumulated in our Milky-Way like galaxy, their buoyancy force overcomes the
magnetic tension of the toroidal disk field. As field lines open up, they
enable anisotropically diffusing CRs to escape into the halo and to accelerate
a bubble-like, CR-dominated outflow. However, these bubbles are invisible in
our simulated gamma-ray maps of hadronic pion-decay and secondary
inverse-Compton emission because of low gas density in the outflows. By
adopting a phenomenological relation between star formation rate (SFR) and
far-infrared emission and assuming that gamma rays mainly originate from
decaying pions, our simulated galaxies can reproduce the observed tight
relation between far-infrared and gamma-ray emission, independent of whether we
account for anisotropic CR diffusion. This demonstrates that uncertainties in
modeling active CR transport processes only play a minor role in predicting
gamma-ray emission from galaxies. We find that in starbursts, most of the CR
energy is "calorimetrically" lost to hadronic interactions. In contrast, the
gamma-ray emission deviates from this calorimetric property at low SFRs due to
adiabatic losses, which cannot be identified in traditional one-zone models.
