How does the supercurrent expression $\vec{j}_s=-\frac{n_se^2}{m}\vec{A}$ arise in Coulomb gauge?

Let's be clear that you are trying to get the London supercurrent ($\vec{j} = - \frac{n_s e^2}{mc}\vec{A}$) which only holds for the wavefunction of a superfluid is rigid and has $ avg(p)=0$ from the GL-theory order parameter.

If you change the phase of the order parameter by $\theta'(r)$ then apply the canonical momentum operator $ P = \frac{\hbar}{i} \nabla + 2 e A$ you get, \begin{equation} P \psi(r) e^{i\theta'(r)} = \Bigg( \frac{\hbar}{i} \nabla + 2 e \bigg(A+\frac{\hbar}{2e} \nabla \theta \bigg) \Bigg)\psi(r) e^{i \theta'(r)}. \end{equation} From this you can see that your choice of gauge transformation is as follows \begin{equation} A(r) \rightarrow A(r) + \frac{\hbar}{2e} \nabla \theta \end{equation} This tells you the vector potential and the phase depend on the choice of the gauge but all other quantities such as free energy and magnetic field are gauge invariant. For a bulk superconductor, the ground state has a constant order parameter which means its magnitude is constant and it's phase only varies very slowly with position r, and this is the so-called phase stiffness. Therefore, you can derive the free energy from the GL theory to be \begin{equation} F_s = F_0 +\rho_s \int d^3r \Bigg( \nabla \theta + \frac{2e}{\hbar}A\Bigg)^2 \end{equation} where $F_0$ is the free energy of the ground state and $\rho_s$ is the superfluid stiffness. Now if we choose the Coulomb gauge $\nabla.A=0$, then there will be a free energy cost if we increase $\nabla \theta$ more. No to minimize the gradient of the free energy you have to take $\theta(r)$ to be constant throughout the superconductor. This is the long-range order in the superconductor. Considering this $\nabla\theta =0$ and the current density becomes \begin{equation} j_s= - \rho_s \bigg(\frac{2e}{\hbar}\bigg)^2 A \end{equation} which is exactly the same as London's current density considering $2\rho_s= |\psi|^2$

I find Feynman's explanation [3rd volume, chapter on Superconductivity] to be very clear. First, he essentially derives the expression $$\vec{j}_s=-\Big[\frac{e\hbar}{m}\nabla\theta+\frac{2e^2}{m}\vec{A}\Big]$$ apart from some factors. Since for the given wavefunction, the probbaility current density $P$ is time independent, he uses the continuity equation to show that $${\rm div}~\vec{j}_s=-\frac{\partial P}{\partial t}=0.$$ Therefore, he obtains $$\nabla^2\theta\propto {\rm div}~\vec{A}.$$ Now, from the vector identity $$\nabla\times\nabla\theta=0$$ and in addition, in the Coulomb gauge $$\nabla^2\theta=\nabla\cdot\nabla\theta=0.$$ Here comes the important bit. There exists no nonzero vector function that is both diveregnce-free and curl-free, and goes to zero at infinity sufficiently rapidly. For the proof, consult Helmholtz’s theorem for vector functions by Peter Young.

Therefore, the only solution is $\nabla\theta=0$ which means that $\theta$ has a spatially uniform profile everywhere! Therefore, it is clear that in the Coulomb gauge $\nabla\theta=0$ and $\vec{j}_s\propto \vec{A}$.

As for me, it seems, that the expression for the current doesn't require choice of gauge. We start from the Hamiltonian for the part of free energy with gradient for Cooper pairs: $$ F = \int \frac{\hbar^2}{4 m} \left|(\nabla - \frac{2 i e}{\hbar} \mathbf{A}) \psi \right|^2 d V $$ The standard procedure for deriving Noether currents, prescribes making a space-dependent transformation $\psi (x) \rightarrow \psi (x) e^{i \alpha (x)}$, $\psi^{*} (x) \rightarrow \psi^{*} (x) e^{-i \alpha (x)}$ $$ \delta F = \int (\nabla \alpha) \frac{e \hbar}{2m} \left(\psi^{*} (\nabla - \frac{2 i e}{\hbar} A) \psi - \psi (\nabla + \frac{2 i e}{\hbar} A) \psi^{*} \right) $$ Which gives simply the aforementioned expression for $j_s$. In this variation, it was unnecessary to impose $\nabla \cdot A = 0$.