If you're using a variational principle this way - which is what Noether's Second Theorem covers - using it with an already-given Lagrangian, then you're seriously under-utilizing the theorem! The second theorem is a model-building technique for deriving action principles and reduces mostly to triviality, once the model is already set up.
The right way to use it is to start out by making no assumptions about what the Lagrangian density is - other than it have the desired symmetry property and then conclude that it has the form that you originally wrote in by hand and started out by assuming, or a generalization thereof.
When used this way with gauge fields, the results are applications of instances of the Utiyama Theorem.
Suppose we have a dynamics given by an action principle with the action integral
$$S = \int 𝔏(q,v,A,V) d^4 x,$$
where the Lagrangian density $𝔏$ depends on fields $\left(q^A: 0 ≤ A < m\right)$ and $\left(A^a_μ: 0 ≤ a < n, 0 ≤ μ < 4\right)$ and their gradients,
$$v^A_μ ≡ ∂_μ q^A,\quad V^a_{μν} ≡ ∂_ν A^a_μ.$$
For simplicity, we assume that $𝔏$ has no explicit dependence on the coordinates $x = \left(x^μ: 0 ≤ μ < 4\right)$. It won't change the analysis in any way, though a more general analysis that also for variations on $x$, too, should also account for variations on $dx$ and work directly with the Lagrangian four-form $L = 𝔏 d^4x$, instead of just with the Lagrangian density $𝔏$.
In here and the following I will use index notation for the components and the Einstein summation convention on indices.
Suppose we have a gauge symmetry
$$δq^C = τ^C_{aB} λ^a q^B,\quad δA^c_μ = -∂_μ λ^c + f^c_{ab} λ^a A^b_μ,$$
that leaves the Lagrangian density invariant. We assume that the coefficients are constant and that their respective matrices
$$τ_a = \left(τ^C_{aB}: 0 ≤ B,C < m\right),\quad f_a = \left(f^c_{ab}: 0 ≤ b,c < n\right),$$
satisfy the matrix relations
$$τ_a τ_b - τ_b τ_a = f^c_{ab} τ_c,\quad f_a f_b - f_b f_a = f^c_{ab} f_c,$$
and that
$$f^c_{ab} = -f^c_{ba}.$$
We also assume the symmetry is transparent with respect to differentiation so that we can write down the corresponding transforms for the gradients:
$$δv^C_μ = ∂_μ\left(δq^C\right) = τ^C_{aB} ∂_μλ^a q^B + τ^C_{aB} λ^a v^B_μ,\\δV^c_{μν} = ∂_ν\left(δA^c_μ\right) = -∂_ν∂_μ λ^c + f^c_{ab} ∂_νλ^a A^b_μ + f^c_{ab} λ^a V^b_{μν}.$$
For convenience, define the following partial derivatives of the Lagrangian density:
$$𝔣_A = \frac{∂𝔏}{∂q^A},\quad 𝔭^μ_A = \frac{∂𝔏}{∂v^A_μ}, \quad 𝔍^μ_a = \frac{∂𝔏}{∂A^a_μ},\quad 𝔊^{μν}_a = \frac{∂𝔏}{∂V^a_{μν}}.$$
Then the variational on the action can be written as:
$$δS = \int δ𝔏 d^4x = \int \left(δq^A 𝔣_A + δv^A_μ 𝔭^μ_A + δA^a_μ 𝔍^μ_a + δV^a_{μν} 𝔊^{μν}_a\right) d^4x,$$
and we assume $δS = 0$, and indeed that $δ𝔏 = 0$.
The invariance is with respect to gauge transforms where the $λ^a$ may be any sufficiently-smooth function of the position. By "sufficiently smooth", we'll assume that they are continuous up to at least the second order of derivatives so that $∂_μ∂_νλ^a = ∂_ν∂_μλ^a$. Then, the integrand in the variational $δS$ must equate to zero - order-by-order in $λ^a$, starting with the highest order derivatives - the second order:
$$0 = δS = \int \left(⋯ + \left(-∂_ν∂_μ λ^a\right) 𝔊^{μν}_a\right) d^4x\quad⇒\quad 0 = ∂_ν∂_μ λ^a 𝔊^{μν}_a = \frac{1}{2} ∂_ν∂_μ λ^a \left(𝔊^{μν}_a + 𝔊^{νμ}_a\right),$$
using the above symmetry condition on the second-order partial derivatives of $λ^a$. By the arbitrariness of $λ^a$, this implies:
$$𝔊^{μν}_a + 𝔊^{νμ}_a = 0\quad⇒\quad 𝔊^{νμ}_a = -𝔊^{μν}_a.$$
Thus, the total differential of $𝔏$ can be written as:
$$d𝔏 = ⋯ + dV^a_{μν} 𝔊^{μν}_a = ⋯ + \frac{1}{2}d\left(V^a_{μν} - V^a_{νμ}\right) 𝔊^{μν}_a,$$
using the just-derived anti-symmetry of $𝔊^{νμ}_a$.
Therefore, $𝔏$ may depend on the gradients $V^a_{μν}$ only through the combination
$$F^c_{μν} ≡ V^c_{νμ} - V^c_{μν} = ∂_μA^c_ν - ∂_νA^c_μ.$$
Thus
$$d𝔏 = ⋯ - \frac{1}{2}dF^c_{μν} 𝔊^{μν}_c.$$
The respective variational in $F^c_{μν}$ is:
$$δF^c_{μν} = δV^c_{νμ} - δV^c_{μν} = f^c_{ab} \left(∂_μλ^a A^b_ν - ∂_νλ^a A^b_μ\right) + f^c_{ab} λ^a F^b_{μν}.$$
Next, for the first order derivatives of $λ^a$, after substituting the variationals $δF^c_{μν}$ for $δV^c_{μν}$, we have:
$$0 = δS = \int \left(⋯ + \left(τ^C_{aB} ∂_μλ^a q^B\right) 𝔭^μ_C + \left(-∂_μ λ^a\right) 𝔍^μ_a - \frac{1}{2} f^c_{ab} \left(∂_μλ^a A^b_ν - ∂_νλ^a A^b_μ\right) 𝔊^{μν}_c\right) d^4x.$$
Using the anti-symmetry of $𝔊^{μν}_c$, we can write the last term as
$$- \frac{1}{2} f^c_{ab} \left(∂_μλ^a A^b_ν - ∂_νλ^a A^b_μ\right) 𝔊^{μν}_c = -f^c_{ab} ∂_μλ^a A^b_ν 𝔊^{μν}_c.$$
Thus,
$$0 = \int ⋯ + ∂_μλ^a \left(τ^C_{aB} q^B 𝔭^μ_C - 𝔍^μ_a - f^c_{ab} A^b_ν 𝔊^{μν}_c\right) d^4x.$$
By the arbitrariness of $λ^a$, the bracketed multiples of $∂_μλ^a$ must vanish, thereby giving us an expression for $𝔍^μ_a$:
$$𝔍^μ_a = τ^C_{aB} q^B 𝔭^μ_C - f^c_{ab} A^b_ν 𝔊^{μν}_c,$$
and reducing the total differential of $𝔏$ to:
$$\begin{align}
d𝔏
&= dq^A 𝔣_A + dv^A_μ 𝔭^μ_A + dA^a_μ 𝔍^μ_a - \frac{1}{2} dF^c_{μν} 𝔊^{μν}_c\\
&= dq^A 𝔣_A + dv^A_μ 𝔭^μ_A + dA^a_μ τ^C_{aB} q^B 𝔭^μ_C - f^c_{ab} dA^a_μ A^b_ν 𝔊^{μν}_c - \frac{1}{2} dF^c_{μν} 𝔊^{μν}_c.
\end{align}$$
Using the anti-symmetry of $f^c_{ab}$ and of $𝔊^{μν}_c$, we can write
$$f^c_{ab} dA^a_μ A^b_ν 𝔊^{μν}_c = \frac{1}{2} f^c_{ab} \left(dA^a_μ A^b_ν + A^a_μ dA^b_ν\right) 𝔊^{μν}_c = \frac{1}{2} d\left(f^c_{ab} A^a_μ A^b_ν\right) 𝔊^{μν}_c.$$
Using integration by parts, and a little index-switching, we can write:
$$dA^a_μ τ^C_{aB} q^B 𝔭^μ_C = d\left(τ^A_{aB} A^a_μ q^B \right) 𝔭^μ_A - dq^A τ^C_{aA} A^a_μ 𝔭^μ_C.$$
Thus,
$$d𝔏 = dq^A \left(𝔣_A - A^a_μ τ^C_{aA} 𝔭^μ_C\right) + d\left(v^A_μ + τ^A_{aB} A^a_μ q^B\right) 𝔭^μ_A - \frac{1}{2} d\left(F^c_{μν} + f^c_{ab} A^a_μ A^b_ν\right) 𝔊^{μν}_c.$$
This shows that the dependence of the Lagrangian density on the field gradients reduces to a dependence on their gauge-covariant derivatives, and that there be no other dependencies on $A^a_μ$. Correspondingly, we rewrite $v$ and $F$ as:
$$
v^A_μ = ∂_μv^A ⇒ v^A_μ = ∂_μq^A + τ^A_{aB} A^a_μ q^B,\\
F^c_{μν} = ∂_μA^c_ν - ∂_νA^c_μ ⇒ F^c_{μν} = ∂_μA^c_ν - ∂_νA^c_μ + f^c_{ab}A^a_μA^b_ν,
$$
and adjust the differential coefficient of $q^A$ accordingly:
$$𝔣_A ⇒ 𝔣_A - A^a_μ τ^C_{aA} 𝔭^μ_C,$$
so that the total differential of $𝔏$ now becomes:
$$d𝔏 = dq^A 𝔣_A + dv^A_μ 𝔭^μ_A - \frac{1}{2} dF^c_{μν} 𝔊^{μν}_c.$$
With these adjustments, all dependency on the first-order derivatives are removed from the variationals and we have:
$$δv^C_μ = τ^C_{aB} λ^a v^B_μ,\quad δF^c_{μν} = f^c_{ab} λ^a F^b_{μν}.$$
For gauge fields $𝔍^μ_a$ is the current density (corresponding to Maxwell's charge density $ρ$ and current density $𝐉$) and $𝔊^{μν}_c$ is the response field (corresponding to Maxwell's $𝐃$ and $𝐇$). The expression derived, above, for $𝔍^μ_a$ shows that there is a gauge-dependent self-interaction involving the response fields and potentials, as well as an extra contribution coming from the $q^A$ field and the response field $𝔭^μ_A$ for the corresponding field strength $v^A_μ$. This gets factored into an adjustment on the expressions for the field strength $v^A_μ$ for $q^A$ and $F^c_{μν}$ (which corresponds to Maxwell's $𝐄$ and $𝐁$) for the potentials $A^a_μ$ (which correspond to the electric/scalar potential $φ$ and magnetic/vector potential $𝐀$), so that the field strengths are now non-linear in terms of the potentials $A^a_μ$ and fields $q^A$.
Finally, at order-zero in $λ^a$, we now have:
$$0 = δS = \int \left(τ^C_{aB} λ^a q^B\right) 𝔣_C + \left(τ^C_{aB} λ^a v^B_μ\right) 𝔭^μ_C - \frac{1}{2} \left(f^c_{ab} λ^a F^b_{μν}\right) 𝔊^{μν}_c d^4x,
$$
which leads to the following identity:
$$τ^C_{aB} q^B 𝔣_C + τ^C_{aB} v^B_μ 𝔭^μ_C = \frac{1}{2} f^c_{ab} F^b_{μν} 𝔊^{μν}_c,$$
which is actually just the continuity equation for the current density.
Edit: The continuity equation, derived here, is off-shell; i.e. independent of any application of the (Euler-Lagrange) field equations. That's a key point of the second Noether Theorem. There's an off-shell continuity equation. How is it to be interpreted? As a pre-condition for the sources in field equation. In particular, the current density associated with the sources (here: the $q^A$ field) must first satisfy the continuity equation before they can quality as sources for the field equations. Otherwise, the field equations are ill-posed.
The Euler-Lagrange equations with the original $𝔣_A$ were:
$$∂_μ 𝔭^μ_A = 𝔣_A,\quad ∂_ν 𝔊^{μν}_a = 𝔍^μ_a.$$
With the modified $𝔣_A$, the first of these becomes:
$$∂_μ 𝔭^μ_A = 𝔣_A + A^a_μ τ^C_{aA} 𝔭^μ_C.$$
The continuity equation for the current density follows by the anti-symmetry in $𝔊^{μν}_a$:
$$∂_μ 𝔍^μ_a = ∂_μ ∂_ν 𝔊^{μν}_a = 0.$$
If you substitute for $𝔍^μ_a$ and use the equations for the modified $𝔣_A$ and for the adjusted $v^a_μ$ and $F^c_{μν}$, you will get the equation just derived from the order-zero variation in $λ^a$.
Further symmetry principles can be used to narrow down on the range of possibilities for Lagrangian density $𝔏$. In particular, requiring $𝔏$ to be locally Lorentz-invariant reduces it to a function of the fields $q^A$ and of the invariant combinations formed from the field gradients, e.g. $g^{μρ}g^{νσ}F^a_{μν}F^b_{ρσ}$, $g^{μν}v^A_μv^B_ν$ and $g^{μρ}g^{νσ}v^A_μF^b_{νρ}v^C_σ$. So, the application of the theorem - which for gauge fields is an instance of the Utiyama Theorem - to all the symmetries, tightly constrains the range of possible Lagrangian densities.
