  We present an in-depth investigation of galaxy clustering based on a new
suite of realistic large-box galaxy-formation simulations in $f(R)$ gravity,
with a subgrid physics model that has been recalibrated to reproduce various
observed stellar and gas properties. We focus on the two-point correlation
functions of the luminous red galaxies (LRGs) and emission line galaxies
(ELGs), which are primary targets of ongoing and future galaxy surveys such as
DESI. One surprising result is that, due to several nontrivial effects of
modified gravity on matter clustering and the galaxy-halo connection, the
clustering signal does not depend monotonically on the fifth-force strength.
For LRGs this complicated behaviour poses a challenge to meaningfully
constraining this model. For ELGs, in contrast, this can be straightforwardly
explained by the time evolution of the fifth force, which means that weaker
$f(R)$ models can display nearly the same -- up to $25\%$ -- deviations from
$\Lambda$CDM as the strongest ones, albeit at lower redshifts. This implies
that even very weak $f(R)$ models can be strongly constrained, unlike with most
other observations. Our results show that galaxy formation acquires a
significant environment dependence in $f(R)$ gravity which, if not properly
accounted for, may lead to biased constraints on the model. This highlights the
essential role of hydrodynamical simulations in future tests of gravity
exploring precision galaxy-clustering data from the likes of DESI and Euclid.
