Difference between covariant derivatives in general relativity and electromagnetism

Your calculation of the gauge covariant derivative is incorrect. The $nabla$ operator changes depending on what it acts upon, just like the covariant derivative in GR. In general a gauge covariant derivative changes based on the representation under which the object it acts upon transforms. This is actually the same as for the covariant derivative in GR as well...there the group is basically $GL(4,mathbb{R})$ and the analog of the vector potential is the Christoffel symbols (to make the two formalisms really precisely identical though one should introduce vierbeins and a spin connection).

If you go and learn about non-abelian gauge theories, this point becomes much clearer because there we speak of fields directly as transforming under some representation or another while in Maxwell theory (which has gauge group $U(1)$), it's the charge that labels the representation so some of this is easily muddled. If you have some familiarity with differential geometry (differentials and so on) and Lie algebras, I can recommend section 2 of this paper (there may be better sources, but this is the first that comes to mind for this particular issue). The goal in that paper is anomalies in quantum field theory, but section 2 is a lightning review of precisely the issues I'm trying to point out here. My recommendation is to try and prove to yourself, using the machinery in the linked paper, how the product rule works out using the fact that a product of fields transforms under the tensor product representation of the representations of the factors in the product (for the GR analog, this is the statement that $a_mu b_nu$ transforms as a tensor). For more introductory sources on the matter, essentially any introductory text on QFT should point these things out in varying levels of detail.

For the sake of this answer, however, let me point out the following: suppose that $a_alpha$ has charge $q_a$, meaning it transforms as $a_alpharightarrow e^{iq_aLambda}a_alpha$ where $Lambda=Lambda(x)$ is the gauge parameter (in the sense $A_murightarrow A_mu+partial_muLambda$ up to conventions about where the gauge coupling appears). Similarly suppose $b_beta$ has charge $q_b$.

It then follows that $a_alpha b_beta$ has charge $q_a+q_b$ since if we were to apply a gauge transformation we would find $$ a_alpha b_betarightarrow e^{iq_aLambda}a_alpha e^{iq_bLambda}b_beta=e^{i(q_a+q_b)Lambda}a_alpha b_beta. $$ It therefore follows that the appropriate covariant derivative acting on the product is (as pointed out in some of the comments under OP, it's a good idea to make explicit where your derivatives act, by the way) $$ D_mu(a_alpha b_beta)=partial_mu(a_alpha b_beta)+i(q_a+q_b)A_mu. $$ The product rule now follows trivially from the product rule of the coordinate derivative and regrouping terms.

Again, I want to stress that this is really exactly the same as the statement that the expression for $nabla_mu V_nu$ and $nabla_mu T_{nulambda}$ are not the same...they have different additions of Christoffel symbols because the index structures are different, and the index structures are identical here to the group representations we are working with.