I think I've figured this out. The point is that, the rigorous meaning one can draw from the formal covariance of $J^\mu$ is that the momentum-space coefficient functions of $J^\mu$ (i.e. the functions in front of monomials of $a_p$ and $a^\dagger_p$) transform covariantly under the change of variable $p\to \Lambda p$. The covariance of the coefficient functions is unaffected by normal ordering, and is sufficient to give rise to the covariance of $:J^\mu:$. The rest of this answer will be an elaboration of the first paragraph.

Let me first clarify the notations used and the meaning of the formal covariance of the ill-defined current $J^\mu$. I'm going to ignore the spin degrees of freedom in this discussion, but one should see the generalization to include spin only involves a straightforward (but perhaps cumbersome) change of notations. I'm also ignoring the spacetime dependence, that is to say I'm only considering the covariance of $J^\mu(0)$, and the generalization to $J^\mu(x)$ is straightforward and easy.

In the context of my question, $U(\Lambda)$ is defined as such that

$$U(\Lambda) a_{p} U^{-1}(\Lambda)=\sqrt{\frac{E_{\Lambda p}}{E_p}}a_{\Lambda p}.$$

The covariance of $J^\mu$ must be understood in a very formal and specific sense, the sense in which the covariance is formally proved. For example, in the case of a fermionic bilinear:

$$U(\Lambda)J^{\mu}U(\Lambda)^{-1}=U\bar{\psi}\gamma^{\mu}\psi U^{-1}\\ =U\bar{\psi}_iU^{-1}(\gamma^{\mu})_{ij}U \psi_j U^{-1}=\bar{\psi}D(\Lambda)\gamma^{\mu}D(\Lambda)^{-1}\psi= \Lambda^{\mu}_{\ \ \nu}\bar{\psi}\gamma^{\nu}\psi, $$

where $D(\Lambda)$ is the spinor representation of Lorentz group, typically constructed via Clifford algebra. Note in this formal proof, what's important is that, under the change $a_{p}\to \sqrt{\frac{E_{\Lambda p}}{E_p}}a_{\Lambda p}$ (ignoring spin indices of course) the elementary field transforms as $\psi \to D(\Lambda)\psi$. In the proof, no manipulation of operator ordering and commutation relations ever occurs: all we do is to do a change of integration variable, and let the algebraic properties of the coefficient functions take care of the rest. In fact, we'd better not mess with the operator ordering, as it can easily spoil the formal covariance (example: $H=\int \text{d}p\frac{1}{2}E_{p}(a_p a_p^\dagger+a_p^\dagger a_p)＝\int \text{d}p E_{p}(a_p^\dagger a_p＋\delta(0))$, see my longest comment under drake's answer).

To explain what's going on in more details without getting tangled with notational nuisances, let me remind you again I'll omit the spin degrees of freedom, but it should be transparent enough by the end of the argument that it's readily generalizable to spinor case, since all that matters is that we know the coefficient functions(even with spin indices) transform covariantly. The mathematical gist is, after multiplying the elementary fields and grouping c/a operators (during the grouping no operator ordering procedure should be performed at all, e.g. $a^\dagger(p_1)a(p_2)$ and $a(p_2)a^\dagger(p_1)$ should be treated as two independent terms), a typical monomial term in $J^\mu(0)$ has the form

$$ \int \left(\prod\limits_{i=1}^{n}\text{d}p_i\right)M(\{a^\dagger(p_i), a(p_i)\})f^\mu(\{p_i\}),$$

where $M$ is a monomial of c/a operators not necessarily normally ordered, but has an ordering directly from the multiplication of elementary fields.

The formal covariance of $J^\mu$ means

$$\Lambda^\mu_{\ \ \nu}\int \left(\prod\limits_{i=1}^{n}\text{d}p_i\right)M(\{a^\dagger(p_i), a(p_i)\})f^\nu(\{p_i\})\\=\int \left(\prod\limits_{i=1}^{n}\text{d}p_i\right)M(\{a^\dagger(\Lambda p_i), a(\Lambda p_i)\})f^\mu(\{p_i\})\\=\int \left(\prod\limits_{i=1}^{n}\text{d}q_i\right)\left(\prod\limits_{i=1}^n \frac{E_{\Lambda^{-1} q_i}}{E_{q_i}}\right) \left(\prod\limits_{i=1}^{m}\sqrt{\frac{E_{q_i}}{E_{\Lambda^{-1} q_i}}}\right) M(\{a^\dagger(q_i), a(q_i)\})f^\mu(\{\Lambda^{-1}q_i\}) ,$$

where $\prod\limits_{i=1}^n {E_{\Lambda^{-1} q_i}}/{E_{q_i}}$ comes from the transformation of measure and $\prod\limits_{i=1}^{m}\sqrt{{E_{q_i}}/{E_{\Lambda^{-1} q_i}}}$ from the transformation of c/a operators in $M$. This is equivalent to

$$f^\mu(\{\Lambda^{-1}q_i\})\left(\prod\limits_{i=1}^n \frac{E_{\Lambda^{-1} q_i}}{E_{q_i}}\right) \left(\prod\limits_{i=1}^{m}\sqrt{\frac{E_{q_i}}{E_{\Lambda^{-1} q_i}}}\right)=\Lambda^\mu_{\ \ \nu}f^\nu(\{q_i\}).$$

The above equation makes completely rigorous sense since it's a statement about c-number functions. Obviously, this equation is sufficient to prove the covariance of the normal ordering

$$ \int \left(\prod\limits_{i=1}^{n}\text{d}p_i\right):M(\{a^\dagger(p_i), a(p_i)\}):f^\mu(\{p_i\}),$$

since on the operator part only a change of integration variable is needed for the proof.

So let's recapitulate the logic of this answer:

1. The current is only covariant when written in a certain way, but not in all ways. (recall the free scalar field Hamiltonian example: $H=\int \text{d}p\frac{1}{2}E_{p}(a_p a_p^\dagger+a_p^\dagger a_p)＝\int \text{d} pE_{p}(a_p^\dagger a_p＋\delta(0))$, which is formally covariant in the first form but not in the second form.)

2. In that certain way where the current is formally covariant, the formal covariance really means a genuine covariance of the coefficient functions.

3. The covariance of the coefficient functions is sufficient to establish the covariance of the normally ordered current.