The most general first-order symmetry operator for the Dirac equation coupled to arbitrary fluxes is derived. Such an operator is given in terms of an inhomogeneous form which is a solution to a coupled system of first-order partial differential equations, the generalized conformal Killing-Yano system.