In a draft answer to another question about gauge transformations, I played around with demonstrating the action of a gauge transformation on the Lagrangian density. Beginning with the classical electromagnetism Lagrangian density \begin{align} \mathcal{L}_{\text{E&M}} &= \frac{\epsilon_0}{2}\mathbf{E}^2 + \frac{1}{2\mu_0}\mathbf{B}^2 - \rho\Phi + \mathbf{J}\cdot\mathbf{A} + \mathcal{L}_{Q}(\rho,\mathbf{J}) \\ & = \frac{\epsilon_0}{2}\left(-\nabla\Phi - \frac{\partial \mathbf{A}}{\partial t}\right)^2 + \frac{1}{2\mu_0}(\nabla\times\mathbf{A})^2 - \rho\Phi + \mathbf{J}\cdot\mathbf{A} + \mathcal{L}_{Q}(\rho,\mathbf{J}), \end{align} where $\mathcal{L}_Q$ is the Lagrangian density for the charge carriers that produce $\rho$ and $\mathbf{J}$. Under a gauge transformation $\Phi\rightarrow\Phi - \partial_t\Lambda$ and $\mathbf{A}\rightarrow \mathbf{A} + \nabla\Lambda$. By construction, the first two terms and $\mathcal{L}_Q$ are gauge invariant [this does not hold in QFT, but let's set that aside for the moment]. Under the gauge transformation, we get \begin{align} \mathcal{L}_{\text{E&M}} & \rightarrow \mathcal{L}_{\text{E&M}} + \rho \partial_t\Lambda + \mathbf{J}\cdot\nabla\Lambda \\ & = \mathcal{L}_{\text{E&M}} + \left[\frac{\partial\rho\Lambda}{\partial t} + \nabla \cdot (\mathbf{J}\Lambda)\right] - \Lambda \left(\frac{\partial \rho}{\partial t} + \nabla\cdot \mathbf{J}\right).\tag1 \end{align} The part of (1) in square brackets is just a collection of surface terms, and the part in large parentheses vanishes by the continuity equation for charge $$\frac{\partial \rho}{\partial t} + \nabla\cdot\mathbf{J} = 0$$ that defines what it means for charge to be locally conserved. So, we can say that as long as charge is locally conserved, the action is invariant under gauge transformations.
Can we not also say, though, that $\Lambda$ plays the role of a Lagrange multiplier that enforces charge conservation? Were I to introduce a Lagrange multiplier to enforce conservation to $\mathcal{L}_{\text{E&M}}$ it would look exactly like (1), up to the surface terms. Admittedly, treating it as such in theories where charge is already conserved is superfluous, but does it produce any interesting results or insights into gauge theory, regardless?
For convenience, the equations of motion that come from the gauge transformed Lagrangian where we treat $\Lambda$ as a Lagrange multiplier are \begin{align} \frac{\delta S}{\delta \Lambda} &=0\Rightarrow & \frac{\partial\rho}{\partial t} + \nabla \cdot \mathbf{J} & = 0 \\ \frac{\delta S}{\delta \Phi} & = 0\Rightarrow & -\epsilon_0\nabla^2\Phi - \nabla \cdot \frac{\partial \mathbf{A}}{\partial t} &= \rho & \Leftrightarrow \nabla\cdot\mathbf{E} & = \frac{\rho}{\epsilon_0} \\ \frac{\delta S}{\delta \mathbf{A}} & = 0\Rightarrow & -\frac{1}{\mu_0} \nabla\times\nabla\times\mathbf{A} - \epsilon_0 \frac{\partial \nabla\Phi}{\partial t} - \epsilon_0 \frac{\partial^2\mathbf{A}}{\partial t^2} & = -\mathbf{J} & \Leftrightarrow \nabla\times\mathbf{B} & = \mu_0\mathbf{J} + \mu_0\epsilon_0 \frac{\partial \mathbf{E}}{\partial t} \\ \frac{\delta S}{\delta \rho} & = 0\Rightarrow & -\Phi + \partial_t\Lambda + \frac{\partial \mathcal{L}_Q}{\partial \rho} &=0 \\ \frac{\delta S}{\delta \mathbf{J}} & = 0 \Rightarrow & \mathbf{A}+ \nabla\Lambda + \frac{\partial \mathcal{L}_Q}{\partial \mathbf{J}} & = 0, \end{align} where derivatives with respect to vector valued quantities are understood to be a form of gradient. Admittedly, those last two equations are not actually equations of motion, since in pretty much every theory of E&M $\rho$ and $\mathbf{J}$ are derived from other variables. Thus, the real equations of motion would either be derived from them using chain rules (with possible space and time derivative shifts), or the relationships have to be inverted to make $\rho$ and $\mathbf{J}$ independent while some other 4-vector field is rendered auxiliary.
