I would like to ask whether fermionic Grassmann fields in a gauge theory path integral (say in QCD) should be chosen to commute or anticommute with ghost and anti-ghost fields. The way most textbooks present it suggests that anti-commutation should be chosen, but I don't know any argument for this. There is also the question of (anti)commutation relations for corresponding operators in the Krein space.
I was asked in the comments where this issue ever came up. It was in derivation of the conserved BRS current in QCD. Let me first explain that in many books it is stated that BRS operator $s$ satisfies graded Leibniz rule with respect to the fermion number, for example $s ( \bar \psi \psi ) = (s \bar \psi ) \psi - \bar \psi s \psi $. On a lecture about anomalies I attended recently it was stated that Leibniz rule graded with respect to the ghost number should be used instead (or at least can be used), so $s ( \bar \psi \psi ) = (s \bar \psi ) \psi + \bar \psi s \psi$ but $s (c^a c^b)=(sc^a) c^b - c^a (s c^b)$. Therefore I am naturally led to consider variations of fields of the form $\phi \to \phi + \epsilon s \phi$, where $\epsilon$ is a Grassmann parameter commuting with $A, \psi, \bar \psi$ but anticommuting with $c^a$ and $\bar c^a$. I found that under these transformations variation of action takes the form $$ \delta S = \int d^4 x (\partial_{\mu} \epsilon) \left[ - F_a^{\mu \nu} D_{\nu} c^a + g \bar \psi \gamma^{\mu} c^a t_a \psi + b^a D^{\mu} c^a - \frac{1}{2} g f_{abc} (\partial^{\mu} \bar c^a) c^b c^c \right]. $$ We see that inside the parenthesis $[]$ we have a conserved current, from now denoted $J_{\mathrm{BRS}}^{\mu}$. After manipulating this current using the equations of motion I found a term $g [\bar \psi, c^a] \gamma^{\mu} t_a \psi$. It turns out that explicit evaluation of the divergence of $J_{\mathrm{BRS}}^{\mu}$ using equations of motion gives zero only if this commutator is taken to be zero. Thus it seems to me that this is the only choice consistent with my choice of the definition of the BRS operator.
Remark I used the Lagrangian $$ \mathcal L = - \frac{1}{4} F^2 + \bar \psi (i \gamma \cdot D - M) \psi + \partial_{\mu} \bar c^a D^{\mu}c^a - A^a_{\mu} \partial^{\mu} b^a + \frac{1}{2} \xi b^2 , $$ with covariant derivative $D=\partial + ig A$.