  We present the first suite of cosmological N-body simulations that
simultaneously include the effects of two different and theoretically
independent extensions of the standard $\Lambda$CDM cosmological scenario -
namely an $f(R)$ theory of Modified Gravity (MG) and a cosmological background
of massive neutrinos - with the aim to investigate their possible observational
degeneracies. We focus on three basic statistics of the large-scale matter
distribution, more specifically the nonlinear matter power spectrum, the halo
mass function, and the halo bias, for which we determine the deviation with
respect to the fiducial $\Lambda$CDM cosmology in the context of both separate
and combined simulations of $f(R)$ MG and massive neutrinos scenarios. Our
results show that while these two extended models separately determine very
prominent and potentially detectable features in all the three statistics, when
we allow them to be simultaneously at work these features are strongly
suppressed, resulting in much weaker deviations from the standard model's
predictions. In particular, when an $f(R)$ gravity model with $f_{R0}=-1\times
10^{-4}$ is combined with a total neutrino mass of $\Sigma _{i}m_{\nu
_{i}}=0.4$ eV, the resulting matter power spectrum, halo mass function, and
bias at z=0 are found to be consistent with the standard model's predictions at
the 10%, 20%, and 5% accuracy levels, respectively. Therefore, our results
imply an intrinsic theoretical limit to the effective discriminating power of
present and future observational data sets with respect to these widely
considered extensions of the standard cosmological scenario in the absence of
independent measurements of the neutrino masses from laboratory experiments,
even though the high-redshift evolution might still allow to partially break
the degeneracy [Abridged].
