Switching to a different formulation can make life easier. To make life easier, I'll do two things:
I'll use non-relativistic quantum field theory (NRQFT), also called second quantization. This allows deriving the general continuity equation as a simple operator equation, without any need for multi-electron wavefunctions or reduced density matrices. The operator equation is derived as a direct consequence of the gauge invariance of the Hamiltonian.
When it's time to introduce a state, I'll use the full state. This is easier than trying to use an explicit reduced density matrix, and for the time-independent continuity equation, the 1P-RDM result falls out automatically anyway. That's because conceptually, when we work with a reduced density matrix, all we're really doing is agreeing to consider only a limited set of operators — such as the set of "one-particle" operators, which includes the electric current operator.
After deriving the desired result, I'll briefly explain how the derivation relates to the von Neumann equation.
The model, part 1: The Hilbert space and Hamiltonian
Let $\varphi(\mathbf{x})$ be a field operator parameterized by $\mathbf{x}$, which denotes the 3-tuple of coordinates for a point in space. Let $\varphi^\dagger(\mathbf{x})$ denote the adjoint of $\varphi(\mathbf{x})$. The OP specifies electrons, which are fermions, so these operators satisfy the anti-commutation relations
$$
\big\{\varphi(\mathbf{x}),\,\varphi^\dagger(\mathbf{y})\big\}
= \delta^3(\mathbf{x} - \mathbf{y})
\tag{1}
$$
and
$$
\big\{\varphi(\mathbf{x}),\,\varphi(\mathbf{y})\big\} = 0
\tag{2}
$$
with $\{A,B\} := AB+BA$. This is QFT's version of saying that the wavefunction must be completely antisymmetric. As suggested in the OP,
I'll pretend that electrons have spin $0$, so the field operators don't need a spin index.
I'll treat the magnetic field as a time-independent external field, expressed in terms of a gauge field $\mathbf{A}(\mathbf{x})$.
For any given choice of the gauge field, the Hamiltonian is
$$
H[\mathbf{A}] =
H_2[\mathbf{A}] +
H_4[\mathbf{A}]
\tag{3}
$$
\begin{align}
H_2[\mathbf{A}]
&= \int d^3x\ \frac{\big(\mathbf{D}\varphi(\mathbf{x})\big)^\dagger
\big(\mathbf{D}\varphi(\mathbf{x})\big)}{2m}
\tag{4}
\\
H_4[\mathbf{A}]
&= \int d^3x\,d^3y\
\big(\varphi(\mathbf{x})\varphi(\mathbf{y})\big)^\dagger
v(\mathbf{x}-\mathbf{y})
\varphi(\mathbf{x})\varphi(\mathbf{y})
\tag{5}
\end{align}
where
$$
\mathbf{D} = \nabla + ie\mathbf{A}(\mathbf{x})
\tag{6}
$$
and where $v(\mathbf{x})$ is a real-valued, decreasing positive function of $|\mathbf{x}|$ that describes the repulsive Coulomb interaction between electrons. For any given gauge field $\mathbf{A}$, the Hamiltonian is manifestly non-negative, and a state $|0\rangle$ that satisfies
$$
\varphi(\mathbf{x})|0\rangle = 0
\tag{7}
$$
for all $\mathbf{x}$ is clearly a state with the lowest possible energy (zero). We can interpret this as the vacuum state, which has no electrons. Applying the operator
$$
\varphi^\dagger(f) := \int d^3x\ f(\mathbf{x})\,\varphi^\dagger(\mathbf{x})
\tag{8}
$$
to any state adds an electron with wavefunction $f$. The Pauli exclusion principle is enforced by the anti-commutation relations, which imply that applying $\varphi^\dagger(f)$ twice to any state gives zero. In words, no two electrons can have the same wavefunction.
The model, part 2: Gauge invariance and observables
We're really defining a whole family of models, one for each choice of the gauge field $\mathbf{A}$. In this context, an "observable" is really a whole family of operators, one operator $F[\mathbf{A}]$ for each choice of the gauge field $\mathbf{A}$, such that the following condition holds for all functions $\theta(\mathbf{x})$:
$$
U[\theta]\,F[\mathbf{A}+\nabla\theta]\,U^{-1}[\theta] = F[\mathbf{A}]
\tag{9}
$$
where the unitary operator $U[\theta]$ is defined by the conditions
\begin{gather}
U[\theta]\,\varphi(\mathbf{x})\,U^{-1}[\theta]
= \exp\big(ie\,\theta(\mathbf{x})\big)\,\varphi(\mathbf{x}).
\tag{10}
\end{gather}
In words, equation (9) says that observables are gauge invariant. An example is
$$
F[\mathbf{A}] = \varphi^\dagger(\mathbf{x})
\left(ie\int_\mathbf{y}^\mathbf{x} \mathbf{A}\cdot d\mathbf{z}
\right)
\varphi(\mathbf{y})
\tag{11}
$$
for a fixed pair of points $\mathbf{x}$ and $\mathbf{y}$. This satisfies the condition (9), so it qualifies as an observable.
Deriving the general (operator) continuity equation
The Hamiltonian $H[\mathbf{A}]$, regarded as a functional of the gauge field, is gauge invariant in the sense defined above. The continuity equation (as an operator equation) can be derived from the gauge invariance of $H[\mathbf{A}]$. To do this, start with the identity
$$
U[\theta]\,H[\mathbf{A}+\nabla\theta]\,U^{-1}[\theta] = H[\mathbf{A}],
\tag{12}
$$
which expresses the gauge invariance of $H[\mathbf{A}]$. Vary both sides with respect to the arbitrary function $\theta$ and then set $\theta=0$ to get
$$
\left(\frac{\delta}{\delta\theta}
H[\mathbf{A}+\nabla\theta]\right)_{\theta=0} +
\left(\frac{\delta}{\delta\theta}
U[\theta]H[\mathbf{A}]U^{-1}[\theta]\right)_{\theta=0} = 0.
\tag{13}
$$
Rewrite the first term using this identity:
$$
\left(\frac{\delta}{\delta\theta}
H[\mathbf{A}+\nabla\theta]\right)_{\theta=0}
\propto \nabla\cdot\mathbf{J}
\tag{14}
$$
with
$$
\mathbf{J}(\mathbf{x})
:= \frac{\delta}{\delta \mathbf{A}(\mathbf{x})}H[\mathbf{A}].
\tag{15}
$$
Equation (15) defines the current-density operator, and by working out the right-hand side, we see that it resembles the expression shown in the OP — except that here it is an operator, involving the electron field operators in place of a single-particle wavefunction. For the second term in (13), use the fact that any operator $X$ satisfies identity
$$
\left(\frac{\delta}{\delta\theta(\mathbf{x})} U[\theta]XU^{-1}[\theta]
\right)_{\theta=0}
\propto [X,J_0(\mathbf{x})]
\tag{16}
$$
where $J_0$ is the charge-density operator
$$
J_0(\mathbf{x}) := e\varphi^\dagger(\mathbf{x})\varphi(\mathbf{x}).
\tag{17}
$$
Altogether, this shows that equation (12), which expresses the gauge invariance of the Hamiltonian, implies the continuity equation
$$
\nabla\cdot\mathbf{J}(\mathbf{x})
\propto [H,J_0(\mathbf{x})]
\tag{18}
$$
for the current- and charge-density operators.
Deriving $\nabla\cdot\mathbf{j}=0$
The continuity equation derived above is an operator equation. It doesn't refer to any state. To answer the question in the OP, we need to introduce a multi-electron state and then consider the "one-particle reduced density matrix." Let $|\psi\rangle$ be an arbitrary $N$-electron state-vector:
$$
|\psi\rangle = \int d^3x_1\cdots d^3x_N\
\Psi(\mathbf{x}_1,...,\mathbf{x}_N)
\varphi^\dagger(\mathbf{x}_1)\cdots
\varphi^\dagger(\mathbf{x}_N)|0\rangle
\tag{19}
$$
where $\Psi$ is the $N$-electron wavefunction. We might as well use the abbreviation
$$
\psi(X) := \frac{\langle\psi|X|\psi\rangle}{\langle\psi|\psi\rangle}
\tag{20}
$$
for any operator $X$, because all predictions that we make with quantum theory ultimately boil down to expressions of this form. I'll follow mathematicians and refer to $\psi(\cdots)$ as the state. It takes any operator as input and returns its expecation value as output. Inserting the operator-valued continuity equation into $\psi(\cdots)$ gives
$$
\nabla\cdot\psi(\mathbf{J})\propto\psi\big([H,J_0]\big).
\tag{21}
$$
If we further assume that $\psi(\cdots)$ is a stationary state — that is, if we assume that $|\psi\rangle$ is an eigenstate of $H$ — then the right-hand side is zero, which leaves
$$
\nabla\cdot\mathbf{j} = 0
\hskip1cm\text{with}\hskip1cm\mathbf{j} := \psi(\mathbf{J}).
\tag{22}
$$
This is the desired result.
Relation to the 1P-RDM
Suppose that $X$ is any operator of the form
$$
X = \int d^3x\,d^3y\ \varphi^\dagger(\mathbf{x})M(\mathbf{x}-\mathbf{y})
\varphi(\mathbf{y})
\tag{23}
$$
where $M$ is an ordinary function or differential operator. The (components of the) current operator $\mathbf{J}$ is one example. For any such operator $X$, the quantity $\psi(X)$ is equal to the trace of $X$ multiplied by the "one-particle reduced density matrix" defined in the OP. This is clear by inspection of equations (19), (20), and (23), taking the anti-commutation relations (1)-(2) into account.
We don't need to make the reduced density matrix explicit, because whenever we talk about a reduced density matrix, all we're really doing is agreeing to consider only a limited set of operators as inputs to the state $\psi(\cdots)$. In the "one-particle" case, we're agreeing to consider only operators of the form (23), which are sometimes called "one-particle" operators because they annihilate-and-then-create just one particle. The operator $\mathbf{J}$ defined above has this form, so the result (22) is automatically a 1P-RDM equation.
Relation to the von Neumann equation
To relate the preceding derivation to the von Neumann equation, start with
$$
\frac{d}{dt}\psi(\cdots) \propto\psi\big([H,\cdots]\big),
\tag{24}
$$
which may be interpreted as the time-evolution equation for the density matrix. (Remember that working with a reduced density matrix really just means that we've agreed to consider only a limited set of operators.) When we write it as in (24), we don't need to worry about how to define an effective one-particle Hamiltonian, which would be awkward because of the two-particle Coulomb interaction. All we need to do is insert the charge-density operator $J_0$ into (24) and use the identity (21) to get
$$
\frac{d}{dt}\psi(J_0) \propto \nabla\cdot\mathbf{j}.
\tag{25}
$$
For a stationary state, the left-hand side is zero, which leaves the desired result (22).
