Things become indeed more complicated in the relativistic case, since both the electron and the proton are spin $1/2$ particles, so the total wave function must be a spinor of $4^2=16$ components. While the spatial part of the wave function can be reduced to an effective one-body form (by the usual introduction of COM and relative coordinates), such a reduction cannot be done exactly for the spinor part. For this reason, the one-body Dirac equation cannot provide exact solutions in case of a finite-mass proton.
This problem was investigated by Salpeter[1] (following the seminal work of Bethe and Salpeter[2]), who derived a wave equation for the hydrogen atom, which was in principle exact on the level of QED (with relativistic and radiative/non-radiative QED effects taken into account). After a series of approximations, resorting to "pure relativistic" corrections (neglecting negative energy, retardation and QED contributions) gives the following two-body wave equation:
$$
\left[\hat{h}_{\text{e}}\otimes I_{\text{P}}+I_{\text{e}}\otimes\hat{h}_\text{P}+V_{\text{CB}}^{++}\right]\Psi(\boldsymbol{r})=E\Psi(\boldsymbol{r}) \ ,
$$
where $\Psi(\boldsymbol{r})$ is a 16-component spinor depending on the relative coordinate $\boldsymbol{r}$,
$$
\hat{h}_{\text{e}}=-i\boldsymbol{\alpha}_\text{e}\cdot\vec{\nabla}+\beta_{\text{e}}m_{\text{e}} \ ,
$$
and
$$
\hat{h}_\text{P}=+i\boldsymbol{\alpha}_\text{P}\cdot\vec{\nabla}+\beta_{\text{P}}m_{\text{P}}
$$
are one-body Dirac Hamiltonians for the electron and proton, respectively, and
$$
V_{\text{CB}}^{++}=\Lambda_{++}\left\{-\frac{Z\alpha}{r}I_{\text{e}}\otimes I_{\text{P}}+\frac{Z\alpha}{2r}\left[\boldsymbol{\alpha}_\text{e}\otimes\boldsymbol{\alpha}_\text{p}+\frac{1}{r^2}(\boldsymbol{\alpha}_\text{e}\cdot\boldsymbol{r})\otimes(\boldsymbol{\alpha}_\text{P}\cdot\boldsymbol{r})\right]\right\}\Lambda_{++}
$$
is the Coulomb-Breit interaction projected onto the positive-energy subspace (also, $\alpha=e^2/(4\pi)$ is the fine structure constant).
Note that Salpeter worked in momentum space; the above formulae are transformed into coordinate space for convenience.
Salpeter treated this equation perturbatively, and obtained finite-mass corrections to the fine structure of hydrogen. The accurate numerical solution of the above equation could also be found with variational techniques, with the small caveat that the $m_{\text{P}}\rightarrow\infty$ limit would not reproduce exactly the one-body Dirac solution. This is due to the presence of the positive energy projection operators (or in the language of QFT, due to the absence of crossed-ladder type Feynman diagrams).
The formula obtained from the perturbative investigation of the Coulomb-Breit interaction reads[3]
$$
E_{nlj}\approx M-\frac{\mu(Z\alpha)^2}{2n^2}-\frac{\mu(Z\alpha)^4}{2n^3}\left[\frac{1}{j+\frac{1}{2}}-\frac{3}{4n}+\frac{1}{4n}\frac{\mu}{M}-(1-\delta_{l0})\left(\frac{1}{j+\frac{1}{2}}-\frac{1}{l+\frac{1}{2}}\right)\frac{\mu^2}{m_{\text{P}}^2}\right] \ .
$$
You can now see that the wikipedia formula is not exact, replacing the electron mass with $\mu$ only works approximately (note e.g. the new $\mu/M$ or $\mu^2/m_{\text{P}}^2$ terms). Such a replacement is more or less justified in "leading order", since the dominant contribution to the fine structure can be calculated by evaluating expectation values of relativistic operators with the non-relativistic wave function (a similar approximation can also be made when calculating finite-mass corrections to the Lamb shift).
The hyperfine structure arising from interactions with the proton spin was left out of $E_{nlj}$ on purpose, even though it is also an ${\cal{O}}(\mu\alpha(Z\alpha)^3\mu/M)$ finite mass contribution. The reason is that treating the proton as a spin $1/2$ elementary particle could only give an incomplete account of hyperfine structure; the naive correction would be calculated with the Dirac $g$-factor $g=2$ instead of the correct $g_{\text{P}}\approx5.586$ value. To correct for this, the excess part of the proton magnetic moment must be introduced via Pauli coupling[4,5].
As a curiosity, note that $E_{nlj}$ depends on the orbital angular momentum quantum number $l$. This means that proton recoil lifts the $j$-degeneracy, and predicts a Lamb-like splitting of the $2s_{1/2}$ and $2p_{1/2}$ levels already on the fine structure level! Of course, this effect is extremely tiny compared to the main QED contribution (for hydrogen, it is ${\cal{O}}(m_{\text{e}}\alpha^4m_{\text{e}}^2/m_{\text{P}}^2)$ against ${\cal{O}}(m_{\text{e}}\alpha^5\ln(\alpha)$).
Update (2023.10.16)
As pointed out by naturallyInconsistent (whom I thank again), there is actually an analitically solvable two-body Dirac Hamiltonian (with Coulomb interaction only) proposed by Marsch[6,7]. After a lengthy derivation, he found (in the rest frame)
$$
{\cal{E}}_{nk}=M\sqrt{1-\frac{2\mu}{M}+\frac{2\mu}{M}f\left(n,k;Z\alpha,\frac{Z\alpha}{2}\right)} \ ,
$$
where
$$
M=m_{\text{e}}+m_{\text{P}} \ \ \ , \ \ \
\mu=\frac{m_{\text{e}}m_{\text{P}}}{m_{\text{e}}+m_{\text{P}}}
$$
are the total and reduced masses, respectively, and
$$
f(p,q;x,y)=\left(1+\frac{x^2}{\left(\sqrt{q^2-y^2}+p-|q|\right)^2}\right)^{-\frac{1}{2}} \ .
$$
We have $k=1$ when the electron and nuclear spin are coupled to a singlet, and either $k=l+1$ (for $J=l+1$) or $k=-l$ (for $J=l-1$) in case of spins coupled to a triplet.
This energy formula describes recoil effects to all orders of $\mu/M$. It is useful to compare the results with those obtained from the "tinkered" one-body Dirac equation of wikipedia:
$$
{\cal{E}}'_{nj}=\mu f\left(n,j+\frac{1}{2},Z\alpha,Z\alpha\right) \ ,
$$
where $m_{\text{e}}$ is manually replaced by $\mu$.
Up to $\mathcal{O}(\alpha^4)$ accuracy, we find
$$
\begin{aligned}
{\cal{E}}_{nk}&\approx M-\frac{\mu(Z\alpha)^2}{2n^2}-\frac{\mu(Z\alpha)^4}{2n^3}\left[\frac{1}{4|k|}-\frac{3}{4n}+\frac{1}{4n}\frac{\mu}{M}\right] \ , \\
{\cal{E}}'_{nj}&\approx \mu-\frac{\mu(Z\alpha)^2}{2n^2}-\frac{\mu(Z\alpha)^4}{2n^3}\left[\frac{1}{j+\frac{1}{2}}-\frac{3}{4n}\right] \ ,
\end{aligned}
$$
to be compared with the perturbative result shown above.
Up to $\mathcal{O}(\mu(Z\alpha)^4\mu/M)$, $E_{nlj}$ agrees perfectly with the formula of Brézin, Itzykson and Zinn-Justin[8]. Apart from the uninteresting rest mass terms, all three formulae produce the same non-relativistic energy (which is not surprising).
The fact that $\mathcal{E}'_{nj}$ misses the higher recoil effects is also expected.
But interestingly, while the $\mathcal{O}(\mu(Z\alpha)^4\mu/M)$ effects match nicely between $\mathcal{E}_{nk}$ and ${E}_{nlj}$, the $\mathcal{O}(\mu(Z\alpha)^4)$ part of $\mathcal{E}_{nk}$ differs from $\mathcal{E}'_{nj}$ and $E_{nlj}$. Due to the effect of nuclear spin, 4$|k|$ seems to take rather different values than $j+1/2$. I will need to think about this. [To be continued...]
References
[1] Salpeter: Phys. Rev. 87 2 328 (1952)
[2] Salpeter, Bethe: Phys. Rev. 84 1232 (1951)
[3] Eides, Grotch, Shelyuto: Phys. Rep. 342 2 63 (2001)
[4] Arnowitt: Phys. Rev. 92 4 1002 (1953)
[5] Newcomb, Salpeter: Phys. Rev. 97 4 1246 (1955)
[3] Marsch: Ann. Phys. 14 5 324 (2005)
[4] Marsch: Ann. Phys. 15 6 434 (2006)
[6] Brézin, Itzykson, Zinn-Justin: Phys. Rev. D 1 2349 (1970)
