Upcoming surveys of cosmic structures will probe scales close to the cosmological horizon, which opens up new opportunities for testing the cosmological concordance model to high accuracy. In particular, constraints on the squeezed bispectrum could rule out the single-field hypothesis during inflation. However, the squeezed bispectrum is also sensitive to dynamical effects of general relativity as well as interactions of matter with residual radiation from the early Universe. In this paper, we present a relativistic simulation pipeline that includes these relativistic effects consistently. We produce light cones and calculate the observed number counts of cold dark matter for five redshift bins between z = 0.55-2.25. We compare the relativistic results against reference Newtonian simulations by means of angular power- and bispectra. We find that the dynamical relativistic effects scale roughly inversely proportional to the multipole in the angular power spectrum, with a maximum amplitude of 10% for ℓ ≲ 5. By using a smoothing method applied to the binned bispectrum we detect the Newtonian bispectrum with very high significance. The purely relativistic part of the matter bispectrum, obtained by subtracting the Newtonian bispectrum from the relativistic one, is detected with a significance of ∼ 3σ, mostly limited by cosmic variance. We find that the pure dynamical relativistic effects accounts for up to 3% and 10% of the total amplitude, respectively in the squeezed and equilateral limits. Our relativistic pipeline for modelling ultra-large scales yields gauge-independent results as we compute observables consistently on the past light cone, while the Newtonian treatment employs approximations that leave some residual gauge dependence. A gauge-invariant approach is required in order to meet the expected level of precision of forthcoming probes of cosmic structures on ultra-large scales.