Gauge transformation on operators Suppose we have two electromagnetic vector potential operators that differ  by a gauge transformation.
$$A'_{\mu}= A_{\mu}-\partial_{\mu}\alpha$$
Now suppose we have
$$\partial_{\mu}\alpha=F(x)\hat a$$
Where $F(x)$ is a function and $\hat a$ an operator.
Is the operator $\hat a$ the identity operator (which means that gauge transformations does not act on the operator level), if not, what is the meaning of the  operator $\hat a$?
 A: There is an example that gauge transformation involves an operator. assume you have a many body wave function $\psi(x_1,..,x_n)$ and you make a gauge transformation to it such that $$\psi(x_1,..,x_n)\rightarrow\exp(-2is\sum_{i<j}\arg(x_i-x_j) )\psi(x_1,..,x_n)$$
in second quantized language this corresponds to $$a^\dagger(x)\rightarrow a^\dagger(x)\exp(-2is\int d^dx \arg(x-x')\hat{\rho}(x')).$$ This can be proved by using the following $$\psi(x_1,...,x_n)=[\prod_i^n\langle vac|a(x_i)][\prod_i^n a^\dagger_{\lambda_i}|vac\rangle]$$ where $\lambda_i$ is the appopriate single particle eigen states. 
If you plug everything in to the second quantized hamiltonian the gauge transformation will read as
$$A\rightarrow A-2s\int d^dx' (\partial_x \arg(x-x'))\hat \rho(x')$$ 
and in functional path integral formalism this will read as
$$A\rightarrow A-2s\int d^dx' \bar\psi(x')(\partial_x \arg(x-x'))\psi(x')$$ where $\psi$ are grassmann numbers
So I showed a specific example that a gauge transformation is in the form of your question.
So I believe in your case a similar situation is happening if you track what is the actual gauge transformation that is made on many body wave function than find what it implies on second quantized operators, finally write your second quantized hamiltonian in terms of those transformed operators you will understand what's happening.
A: The action of classical electromagnetic field is $$\mathcal{L}=-\frac{1}{4}F_{\mu\nu}F^{\mu\nu}=\frac{1}{2}(\|\mathbf{E}\|^{2}-\|\mathbf{B}\|^{2}).$$
The canonical momentum of $A^{i}$ is $$\Pi^{i}=\frac{\partial\mathcal{L}}{\partial\dot{A}_{i}}=-\dot{A}^{i}+\partial^{i}A^{0}=E^{i}.$$
Notice that $A^{0}$ appears in the Lagrangian density of a Lagrange multipier who has no dynamics.
Before fixing gauge redundancies, one can define the Poisson bracket $$\left\{F,G\right\}_{PB}=\sum_{i=1}^{3}\int d^{3}x\left\{\frac{\delta F}{\delta A^{i}(\mathbf{x})}\frac{\delta G}{\delta\Pi^{i}(\mathbf{x})}-\frac{\delta F}{\delta\Pi^{i}(\mathbf{x})}\frac{\delta G}{\delta A^{i}(\mathbf{x})}\right\},$$
on the naive phase space of photons. So the fundamntal Poisson bracket of electromagnetic fields reads $$\left\{A_{i}(\mathbf{x}),\Pi_{j}(\mathbf{y})\right\}_{PB}=\delta_{ij}\delta(\mathbf{x}-\mathbf{y}).$$
Then, the Hamiltonian is
\begin{align}
H_{\mathrm{em}}&=\int d^{3}x\left(\Pi^{i}\dot{A}_{i}-\mathcal{L}_{\mathrm{em}}\right) \\
&=\int d^{3}x\left(\Pi_{i}(\partial^{0}A^{i}-\partial^{i}A^{0})+\Pi_{i}\partial^{i}A^{0}-\mathcal{L}_{\mathrm{em}}\right) \\
&=\int d^{3}x\left(-\Pi_{i}E^{i}+\partial_{i}(\Pi^{i}A^{0})-A^{0}\partial_{i}\Pi^{i}-\mathcal{L}_{\mathrm{em}}\right) \\
&=\int d^{3}x\left(\frac{1}{2}\|\mathbf{E}\|^{2}+\frac{1}{2}\|\mathbf{B}\|^{2}-A^{0}\nabla\cdot\mathbf{E}\right) \\
&=\int d^{3}x\left(\frac{1}{2}\|\mathbf{\Pi}\|^{2}+\frac{1}{2}\|\nabla\times\mathbf{A}\|^{2}-A^{0}\nabla\cdot\mathbf{\Pi}\right),
\end{align}
where $A_{0}$ appears as an auxiliary field which produces the Gauss constraint $\nabla\cdot\mathbf{E}=0$.
The gauge redundancy you are talking about is generated by the functional $$\mathcal{G}[\lambda]=\int d^{3}x\,\alpha(\mathbf{x})\nabla\cdot\mathbf{\Pi}(\mathbf{x}),$$
because
\begin{align}
&\delta_{\mathcal{G}}A_{i}(\mathbf{y})=\left\{A_{i}(\mathbf{y}),\mathcal{G}\right\}_{PB}=\partial_{i}\alpha(\mathbf{y}), \\
&\delta_{\mathcal{G}}\Pi_{i}(\mathbf{y})=\left\{\Pi_{i}(\mathbf{y}),\mathcal{G}\right\}_{PB}=0. 
\end{align}
Indeed, from the above equations, one finds that the functional $\mathcal{G}[\alpha]$ generates an infinitesimal gauge transformation.
Canonical Quantization:
The naive canonical quantization attained is by replacing the classical Poisson bracket by a canonical commutation relation $$\left[\hat{A}_{i}(\mathbf{x}),\hat{\Pi}_{j}(\mathbf{y})\right]=-i\hbar\delta_{ij}\delta(\mathbf{x}-\mathbf{y}).$$
Notice that this naive canonical commutation relation is incompatible with the Gauss constraint $\nabla\cdot\hat{\mathbf{E}}=0$. This is precisely because one has not yet fixed the gauge. i.e the phase space is too large. Indeed, one can use the above commutation relation and the Baker–Campbell–Hausdorff formula and check that
$$e^{i\hat{\mathcal{G}}[\alpha]}\hat{A}_{i}(\mathbf{x})e^{-i\hat{\mathcal{G}}[
\alpha]}=\hat{A}_{i}(\mathbf{x})+\partial_{i}\alpha(\mathbf{x}),$$
which is exactly the operator equation you are talking about.
The $U(1)$-valued operator $$e^{i\hat{\mathcal{G}}[\alpha]}=\exp\left(i\int d^{3}x\,\alpha(\mathbf{x})\nabla\cdot\hat{\mathbf{\Pi}}(\mathbf{x})\right)$$
generates finite gauge transformations on $\hat{\mathbf{A}}(\mathbf{x})$. Since physical states should be gauge invariant, one has to impose the condition $$e^{i\hat{\mathcal{G}}[\alpha]}|\Psi\rangle=|\Psi\rangle$$
for any physical state $|\Psi\rangle$.
