  We investigate the properties of dark matter haloes and subhaloes in an
$f(R)$ gravity model with $|f_{R0}|=10^{-6}$, using a very high-resolution
N-body simulation. The model is a borderline between being cosmologically
interesting and yet still consistent with current data. We find that the halo
mass function in this model has a maximum 20% enhancement compared with the
$\Lambda$CDM predictions between $z=1$ and $z=0$. Because of the chameleon
mechanism which screens the deviation from standard gravity in dense
environments, haloes more massive than $10^{13}h^{-1}M_\odot$ in this $f(R)$
model have very similar properties to haloes of similar mass in $\Lambda$CDM,
while less massive haloes, such as that of the Milky Way, can have steeper
inner density profiles and higher velocity dispersions due to their weaker
screening. The halo concentration is remarkably enhanced for low-mass haloes in
this model due to a deepening of the total gravitational potential. Contrary to
the naive expectation, the halo formation time $z_f$ is later for low-mass
haloes in this model, a consequence of these haloes growing faster than their
counterparts in $\Lambda$CDM at late times and the definition of $z_f$.
Subhaloes, especially those less massive than $10^{11}h^{-1}M_\odot$, are
substantially more abundant in this $f(R)$ model for host haloes less massive
than $10^{13}h^{-1}M_\odot$. We discuss the implications of these results for
the Milky Way satellite abundance problem. Although the overall halo and
subhalo properties in this borderline $f(R)$ model are close to their
$\Lambda$CDM predictions, our results suggest that studies of the Local Group
and astrophysical systems, aided by high-resolution simulations, can be
valuable for further tests of it.
