  We have carried out a set of cosmological hydrodynamical simulations that
follow galaxy formation in $f(R)$ modified gravity models. Our simulations
employ the Illustris-TNG full physics model and a new modified gravity solver
in the AREPO code. For the first time we are able to investigate the degeneracy
in the matter power spectrum between the effects of $f(R)$-gravity and feedback
from active galactic nuclei (AGN), and the imprint of modified gravity on the
properties of galaxies and on the distribution of dark matter, gas and stars in
the universe. $f(R)$-gravity has an observable effect on the neutral hydrogen
power spectrum at high redshift at a level of 20%. For both the F6 and F5
models, this is significantly larger than the predicted errors for the SKA1-MID
survey, making this probe a powerful test of gravity on large scales. A similar
effect is present in the power spectrum of the stars at high redshift. We also
show that rotationally supported disc galaxies can form in $f(R)$-gravity, even
in the partially screened regime. Our simulations indicate that there might be
more disc galaxies in F6 compared to GR, and fewer in F5. Finally, we show that
the back reaction between AGN feedback and modified gravity in the matter power
spectrum is not important in the F6 model but has a sizeable effect in F5.
