Conservation law of colour current in Yang-Mills theories

I think there's a bit of confusion in how you are writing the equation of motion.

TL;DR: the current usually used for the conservation of colour currents is the quark current, i.e. the one referring to the matter component of the lagrangian (as opposed to the gluon gauge field part). This current is not the same current that you'd get from Noether's theorem. And technically Noether's theorem only applies to global symmetries which is not the case for QCD.

---

The equation of motion for the gluon field $F^a_{munu}$ is: $$tag{1} partial^mu F^a_{munu}(x) + f_{abc}A^mu_bF^c_{munu}(x) = - color{red}{j}^a_nu(x),$$ where the lower case $j$ is used for the matter currents, in this case the colour currents of the quarks: $$ j^a_nu(x)= barpsi(x)gamma_nu T_a psi(x) = bar psigamma_nu frac{lambda_a}{2}psi, $$ where $T^a$ are the generators of $SU(3)$ and $lambda_a$ the Gell-Mann matrices.

Now.

In eq. 1, bring the $f_{abc}...$ bit on the RHS and you get: $$tag{2} partial^mu F^a_{munu}(x) = - f_{abc}A^mu_bF^c_{munu}(x) - color{black}{j}^a_nu(x) = color{red}{J}^a_nu(x).$$

Now this $J^a_mu = - f_{abc}A^mu_bF^c_{munu}(x) - color{black}{j}^a_nu(x)$ is:

• The current that appears in the differential form: $$ partial^mu F^a_{munu}(x) = color{black}{J}^a_nu(x) quad leftrightarrow quad mathrm{d}F = J$$
• This current is the "Noether" current. Noether's (first) theorem only applies to global symmetries, whereas QCD is a local $SU(3)$ symmetry so Noether's formalism wouldn't strictly speaking apply as strongly.
  But if you assumed a Yang-Mills lagrangian $$ mathcal{L}_{text{YM}} = mathcal{L}_{text{field}} + mathcal{L}_{text{matter}}$$ and apply the usual formula for the Noether current $$ J^mu = frac{delta mathcal{L}}{delta(partial_mu varphi_i)}delta varphi_i,$$ you'd get: $$ J^mu propto delta mathcal{L}_{text{YM}} propto deltamathcal{L}_{text{field}} + deltamathcal{L}_{text{matter}},$$ i.e. two things, that correspond to the two bits in $J^a_mu$ above the bullet points.
  
  And, nicely, we confirm that the current associated with the matter part of the Yang-Mills Lagrangian is indeed $j^a_mu$ as we had mentioned at the beginning of the answer.
  
  Coming from Noether's theorem, this $J^mu$ is also conserved according to: $$ partial_mu J^mu = 0.$$

---

So, going back to the matter current $j^a_nu$. Is it "covariantly" conserved?

Luckily, we can start form eq. 1 and use the covariant derivative: $$ D^{ab}_mu = delta^{ab}partial_mu + f_{abc}A^{c}_mu $$ to re-write eq. 1 as: $$ D^mu F^a_{munu} = -j_nu^a(x),$$ so the same thing as your third equation but with lower case $j$ i.e. the matter current (quarks).

And, as you have yourself showed, you do end up with: $$ D_mu j^mu_a =0,$$ so yes, the matter current is "covariantly" conserved.

But now you may say "what if I write $D^mu$ as $partial^mu + dots$, where then $partial^mu j^a_mu =0$ and I am left with the other bit".

The justification for $partial^mu j^a_mu =0$ would be another Noether's theorem but only applying to the matter part of the lagrangian. So if you only consider that bit, then sure just be content with $partial^mu j^a_mu =0$ giving you conservation of colour currents.

But if you want to bring in the covariant derivative, then you also have to consider the gauge field part of the lagrangian and then consider the $J^mu$ "full" current discussed above.

---

And as per your GR connection at the very end, note that GR is not a Yang-Mills theory so you cannot as easily draw parallels between the two. See the end of this answer for a more quantitative discussion about this point though.