I am trying to understand the structure of the fermions in non-abelian gauge theories. Disclaimer: my question might be very trivial (I suspect the answer could simply be "a change of basis"), but I would be grateful is someone could shed some light if there's something deeper lurking around the corner.
Let's consider the Yang-Mills Lagrangian
$$ \mathcal{L} = -\frac{1}{4} (F_{\mu \nu}^a)^2 + \bar{\psi}(i \gamma^\mu D_\mu -m)\psi $$
where $ D_\mu = \mathbf{1} \partial_\mu - i g A_\mu^a t^a$, then $\psi$ needs to have both Dirac ($\mu$) and colour (a) degrees of freedom. I get confused when I change the regular approach to the problem and I am not sure if it is a real issue or if I'm just overcomplicating/overthinking this.
For instance, let's consider $SU(N)$ YM where the generators $t^a$ are $N \times N$ matrices, and there are $N^2-1$ of them. Therefore, the covariant derivative $D_\mu$ is an $N \times N$ matrix, and the index $a$ above runs over $N^2-1$ value.
When we start by contracting $D_\mu$ with the $\gamma$-matrices,
$$ (\gamma^\mu D_\mu)_{ij} = \delta_{ij} \gamma^\mu \partial_\mu - i g (\gamma^\mu A_\mu^a) (t^a)_{ij} $$
one gets an $N \times N$ matrix of $4 \times 4$ matrices. The corresponding $4N$-component object this matrix acts upon is the $N$ Dirac spinors arranged in a column.
However, notice that if we start in the following way instead:
$$ (D_\mu \psi)_i = \partial_\mu \psi_i - i g A_\mu^a (t^a)_{ij} \psi_j $$
we get that the covariant derivative acts on $\Psi \equiv (\psi_1, \cdots, \psi_N)$. For all we know, $\Psi$ has no spinor structure since we haven't contracted with the gamma matrices yet.
Contracting with the $\gamma^\mu$, one gets a $4 \times 4$ matrix of $N \times N$ matrices that encodes the same information as before. This time, it would seem that we only have one Dirac spinor, where each component is a $N$-valued singlet.
However, it seems that the matrix $\gamma^\mu D_\mu$ and the $4N$-component spinor look different, even though all we did was to change the order in which we constructed things.
In the first case, we get $N$ spinors corresponding to the $N$ colors of the adjoint representation. In $SU(3)$, this would be like saying we effectively have 3 spinors which correspond to the red/blue/green colors as $(\psi_R, \psi_B, \psi_G)$. In the second case, this identification fails since we only have one big complicated object.
What went wrong? Is this difference simply a change of basis for $\psi$? Is there something relevant we can learn from looking at the YM Lagrangian in these two different ways?
Also, in the first case, when we get $N$ spinors, I am confused about their signification. I always assumed that in QCD, $\psi$ would correspond to a quark, which is a fermion by itself. Does it mean that quarks are fermions that can be described by fermionic fields/degrees of freedom that we call color?