  We compute the QCD$\times$QED (${\cal{O}}(\alpha_s \alpha)$) mixed and
QED$^2$ (${\cal{O}}(\alpha^2)$) corrections to the production of an on-shell
$Z$ boson in hadronic collisions. We obtain them by profiting from the
calculation of the pure QCD terms after taking the corresponding abelian
limits. Therefore, we extend the available knowledge up to complete
next-to-next-to leading order precision in QCD$\oplus$QED. We present explicit
results for the perturbative coefficients and perform the phenomenological
analysis at different collider energies with particular emphasis on the mixed
corrections. We study the contribution from the different channels and discuss
the scale dependence stabilisation effect. We consider a factorisation
approximation for the mixed order terms and show that it fails to reproduce the
exact result. We find that the contributions are small, typically at the few
per mille level, but that under some kinematical conditions they can compete
with the pure QCD NNLO ones.
