As already mentioned in the comment, the exact version was proven in the link. However, if you insisted to show it infinitesimally in your notation, here is how.
First, consider the variation (warning: sloppy of notation use in the following but it is clear in the context I guess)
\begin{equation} \label{delL}
\delta L = -\frac{1}{2} \sum_a \left( \partial_\mu A_\nu^a - \partial_\nu A^a_\mu + gf^{abc}A^b_\mu A^c_\nu \right) \left( \partial_\mu \delta A^a_\nu + gf^{abc}A^b_\mu \delta A^c_\nu - (\mu \leftrightarrow \nu) \right),
\end{equation}
where
\begin{equation} \label{delA}
\delta A^a_\mu = \frac{1}{g} \partial_\mu A^a - f^{abc} \alpha^b A^c_\mu.
\end{equation}
Now, if you substitute \tag{delA} into \eqref{delL} and patiently listed every terms, you will encounter the terms you noticed already:
\begin{equation}
\delta L = -\frac{1}{2}\sum_a F_{\mu\nu}^a
\left[ \partial_\mu \left( \frac{1}{g}\partial_\nu \alpha^a - f^{abc}\alpha^b A^c_\nu \right) + f^{abc}A^b_\mu \left( \frac{1}{g}\partial_\nu\alpha^c - f^{cde}\alpha^d A^e_\nu \right) - (\mu\leftrightarrow\nu) \right]
\end{equation}
I will not show in detail every step, but describe the computation. Expanding the first big round bracket inside the big square bracket, you will have the derivative hitting on the (infinitesimal) gauge parameter,
$
\frac{1}{g}\partial_\mu \partial_\nu \alpha^a - (\mu\leftrightarrow\nu).
$
The above two terms (notice the term with $\mu,\nu$ exchanged actually cancel each other, because
\begin{equation}
\partial_\mu \partial_\nu \alpha^a = \partial_\nu \partial_\mu \alpha^a,
\end{equation}
i.e. partial derivatives commute.
There are also terms,
$
-f^{abc} (\partial_\mu \alpha^b) A^c_\nu + f^{abc}A^b_\mu \partial_\nu \alpha^c
- (\mu\leftrightarrow\nu),
$
where the above first term comes from the first round bracket and the above second term from the second round bracket inside the big square bracket. One can easily show that they cancel each other after renaming the indices and considering the anti-symmetric properties of the spacetime indices $\mu,\nu$ and the internal indices $a,b,c$.
Finally, the most tricky part is the remaining terms,
\begin{equation}
-f^{abc} \alpha^b \partial_\mu A^c_\nu - f^{abc}f^{cde} \alpha^d A^b_\mu A^e_\nu
- (\mu \leftrightarrow \nu).
\end{equation}
The second term above together with their $- (\mu \leftrightarrow \nu)$ counter part can be rewritten using Jacobi identity and then combine with the rest term to get
\begin{equation}
f^{abc}\alpha^b F_{\mu\nu}^c
\end{equation}
Substitute this back to the original variation equation \eqref{delL}, you will get
\begin{equation}
f^{abc} \alpha^b F_{\mu\nu}^c F^{\mu\nu\, a}.
\end{equation}
You may want to verify that the last pieces in this form vanish (relatively) manifest.
