The variational method, within the Hamiltonian formalism of reformulated QED, is used to determine relativistic wave equations for systems of three fermions of arbitrary mass interacting electromagnetically. The derived interaction kernels of the equations are, in essence, the invariant matrices in the lowest order. The equations are used to obtain relativistic O(α2) corrections to the non-relativistic ground-state energy levels of the positronium and muonium negative ions (e+e−e−, μ+e−e−), as well as H−, using approximate variational three-body wavefunctions. The results are compared with other calculations, where available.