  We quantify the impact of galaxy formation on dark matter halo shapes using
cosmological simulations at redshift $z=0$. The haloes are drawn from the
IllustrisTNG project, a suite of magneto-hydrodynamic simulations of galaxies.
We focus on haloes of mass $10^{10-14} M_\odot$ from the 50-Mpc (TNG50) and
100-Mpc (TNG100) boxes, and compare them to dark matter-only (DMO) analogues
and other simulations e.g. NIHAO and Eagle. We further quantify the prediction
uncertainty by varying the baryonic feedback models in a series of smaller 25
Mpc $h^{-1}$ boxes. We find that: (i) galaxy formation results in rounder
haloes compared to the DMO simulations, in qualitative agreement with past
hydrodynamic models. Haloes of mass $\approx 2\times 10^{12} M_\odot$ are most
spherical, with an average minor-to-major axis ratio of $\left< s \right>
\approx 0.75$ in the inner halo, an increase of 40 per cent compared to their
DMO counterparts. No significant change in halo shape is found for low-mass
$10^{10} M_\odot$ haloes; (ii) stronger feedback, e.g. increasing galactic wind
speed, reduces the impact of baryons; (iii) the inner halo shape correlates
with the stellar mass fraction, which can explain the dependence of halo shapes
on different feedback models; (iv) the fiducial and weaker feedback models are
most consistent with observational estimates of the Milky Way halo shape. Yet,
at fixed halo mass, very diverse and possibly unrealistic feedback models all
predict inner halo shapes that are closer to one another than to the DMO
results. This implies that a larger observational sample would be required to
statistically distinguish between different baryonic prescriptions due to large
halo-to-halo variation in halo shapes.
