On the track you're on, $⟨AB⟩$ can indeed differ from $⟨A⟩⟨B⟩$ even if $[A,B]=0$. For a trivial counterexample, use $A=B=x$ on $\mathcal H=L_2(\mathbb R)$, for which in general
$$
⟨AB⟩=⟨x^2⟩>⟨x⟩^2=⟨A⟩⟨B⟩
$$
unless the state is a delta function at the origin.
On the other hand, one expects any integral of motion to be conserved on both real and imaginary time evolution. This follows from the fact that
$$
n(t)=⟨\psi|e^{iHt}\hat Ne^{-iHt}|\psi⟩
$$
should be an analytical function around the real axis, and if it is constant on the real axis it must be constant everywhere.
More generally, if $[\hat N,H]=0$ then you expect any evolution under $H$, under both real and imaginary time, to preserve the eigenstates of $N$. Thus, if $\hat N|\psi⟩=N|\psi⟩$ then
$$
\hat Ne^{iH\tau}|\psi⟩=Ne^{iH\tau}|\psi⟩
$$
for all $\tau\in\mathbb C$. The time evolution under $\tau$ will then also preserve the expectation value of $\hat N$.
There is a problem, on the other hand, if you want to look at the expectation value of $\hat N$ in a superposition of shared eigenstates of $\hat N$ and $\hat H$. To see why, take $|\phi_1⟩$ and $|\phi_2⟩$ to be such shared eigenstates with $\hat N|\phi_i⟩=N_i|\phi_i⟩$ and $\hat H|\phi_i⟩=E_i|\phi_i⟩$, with $N_1\neq N_2$ and $E_1<E_2$. Then, if you propagate along the positive imaginary time axis, you get
\begin{align}
|\psi(it)⟩
&=e^{-iH(it)}\frac1{\sqrt{2}}\left(|\phi_1⟩+|\phi_2⟩\right)
\\&=\frac{e^{E_1t}}{\sqrt{2}}|\phi_1⟩+\frac{e^{E_2t}}{\sqrt{2}}|\phi_2⟩
\end{align}
and therefore
\begin{align}
⟨\psi(it)|\hat N|\psi(it)⟩
&=\frac{e^{2E_1t}}{2}N_1+\frac{e^{2E_2t}}{{2}}N_2
\end{align}
so that
\begin{align}
⟨\hat N⟩:=
\frac{⟨\psi(it)|\hat N|\psi(it)⟩}{⟨\psi(it)|\psi(it)⟩}
=\frac{
e^{2E_1t}N_1+e^{2E_2t}N_2
}{
e^{2E_1t}+e^{2E_2t}
}
=\frac{
N_1+e^{-2(E_1-E_2)t}N_2
}{
1+e^{-2(E_1-E_2)t}
}.
\end{align}
In particular $⟨\hat N⟩$ changes from $(N_1+N_2)/2$ at time $t=0$ to $N_1$ in the limit $it\to i\infty$. (Note, here, that the heuristic above is not violated, since $⟨\hat N⟩(t)=⟨\psi(t)|\hat N|\psi(t)⟩$ is oscillatory - and not constant - for real times.)
This is sort of your case, but the differences only make the non-conservation useful for you.
we want to find a ground state of a system with a hamiltonian commuting with total number of particles.
The above guarantees that (i) the ground state produced by the procedure will be an eigenstate of $\hat N$, (ii) if you specify at the beginning a particle number that you want in your ground state (i.e. by choosing your initial guess in that eigenspace, which is easy to do (and in fact it is very hard not to do)) then the ground state produced by the procedure will retain that eigenstate, and (iii) if for some reason you mess up and you specify a seed with components in more than one $\hat N$ eigenspace, the imaginary-time procedure will automatically select for you the correct ground state, $\hat N$ eigenstate and all.
