In the case of Abelian symmetry, the covariant derivative is defined as $D_\mu\equiv \partial_\mu + ieA_\mu$, where $e$ is an arbitrary constant and the vector field, $A_\mu$ is a called a gauge field. The commutator is $$\left[D_\mu, D_\nu\right]\phi=D_\mu D_\nu\phi-D_\nu D_\mu\phi$$ $$=\left(\partial_\mu + ieA_\mu\right)\left(\partial_\nu + ieA_\nu\right)\phi-\left(\partial_\nu + ieA_\nu\right)\left(\partial_\mu + ieA_\mu\right)\phi$$ $$=\partial_\mu\partial_\nu\phi+ieA_\mu\partial_\nu\phi +ie\left(\partial_\mu A_\nu\right)\phi+ie A_\nu\partial_\mu\phi-\color{red}{e^2A_\mu A_\nu\phi}\tag{a}$$ $$-\partial_\nu\partial_\mu\phi-ieA_\nu\partial_\mu\phi-ie\left(\partial_\nu A_\mu\right)\phi-ieA_\mu\partial_\nu\phi+\color{red}{e^2A_\nu A_\mu\phi}\tag{b}$$ $$=ie\left(\partial_\mu A_\nu-\partial_\nu A_\mu\right)\phi=ieF_{\mu\nu}\phi$$ Where $F_{\mu\nu}=\partial_\mu A_\nu-\partial_\nu A_\mu$. In going from $(\mathrm{a})$ to $(\mathrm{b})$ I can see that 6 out of the 10 terms cancel but why should the terms marked red cancel?
Put another way, I think that $$\left[D_\mu, D_\nu\right]\phi=ieF_{\mu\nu}\phi-e^2\left[A_\mu,A_\nu\right]\phi\tag{c}$$ Now the default answer to this may be 'it is because it is Abelian', but the Abelian nature of the (internal) symmetry transformation (leaving the Lagrangian invariant), $\phi(x)\to e^{i\theta}\phi(x)$ is such that any pair of symmetry transformations commute, $e^{i\theta_1}e^{i\theta_2}=e^{i\theta_2}e^{i\theta_1}$, where $\theta$ is an arbitrary real number. But in the text I'm reading I see no mention of the gauge field $A_\mu$, being commutative.
Now for the non-Abelian case, the situation is very similar, except the covariant derivative is defined as $D_\mu\equiv \partial_\mu + igA_\mu$, where $g$ is an arbitrary constant. As can be seen below, the commutator of the gauge fields is non zero: $$\left[D_\mu, D_\nu\right]\phi=D_\mu D_\nu\phi-D_\nu D_\mu\phi$$ $$=\left(\partial_\mu + igA_\mu\right)\left(\partial_\nu + igA_\nu\right)\phi-\left(\partial_\nu + igA_\nu\right)\left(\partial_\mu + igA_\mu\right)\phi$$ $$=\partial_\mu\partial_\nu\phi+igA_\mu\partial_\nu\phi +ig\left(\partial_\mu A_\nu\right)\phi+ig A_\nu\partial_\mu\phi-\color{blue}{g^2A_\mu A_\nu\phi}\tag{d}$$ $$-\partial_\nu\partial_\mu\phi-igA_\nu\partial_\mu\phi-ig\left(\partial_\nu A_\mu\right)\phi-igA_\mu\partial_\nu\phi+\color{blue}{g^2A_\nu A_\mu\phi}\tag{e}$$ $$=ig\left(\partial_\mu A_\nu-\partial_\nu A_\mu\right)\phi-g^2\left[A_\mu,A_\nu\right]\phi$$
In short, what part of the symmetry transformation for the field, $\phi(x)\to e^{i\theta}\phi(x)$ suggests that the gauge field in eqn. $(\mathrm{c})$, $A_\mu$ is commutative?
So the gauge is such that $A_\mu\to A_\mu-\partial_\mu\lambda$ where $\lambda$ is an arbitrary scalar function of spacetime, but how does this gauge transformation imply that $\left[A_\mu,A_\nu\right]=0$ for the Abelian case?
I have typed details from the notes in this previous question that pertains exactly to the gauge invariance that is assumed in this question.
These sources of information were obtained from notes on QFT from ICL dept. of physics. If anyone needs more information please let me know, thanks!