The inverse of the Faddeev-Popov operator plays a pivotal role within the Gribov-Zwanziger approach to the quantization of Euclidean Yang-Mills theories in Landau gauge. Following a recent proposal (Capri etal., 2014), we show that the inverse of the Faddeev-Popov operator can be consistently coupled to quark fields. Such a coupling gives rise to a local action while reproducing the behaviour of the quark propagator observed in lattice numerical simulations in the non-perturbative infrared region. By using the algebraic renormalization framework, we prove that the aforementioned local action is multiplicatively renormalizable to all orders.