One method is to look at the interaction energy between two charges sources. This can be derived from the path integral
\begin{equation}
Z[J] = e^{i W[J]} = \int DA_\mu e^{i \int d^4 x \left(-\frac{1}{4} F^{\mu\nu}F_{\mu\nu} + A_\mu J^\mu\right)}
\end{equation}
The integral can be done exactly by completing the square because it is a Gaussian integral and yields
\begin{equation}
W[J] = -\frac{1}{2} \int d^4 x \int d^4 y J^\mu(x) G_{\mu\nu}(x, y) J^\nu(y)
\end{equation}
where the propagator is
\begin{equation}
G_{\mu\nu} = \frac{\eta_{\mu\nu}}{\square - i \epsilon} = \int \frac{d^4 p}{(2\pi)^4}\frac{\eta_{\mu\nu} e^{i p_\mu (x^\mu-y^\mu)}}{p^2 - i \epsilon}
\end{equation}
There is a subtlety here because of gauge invariance. Strictly speaking, one cannot invert the differential operator $O_{\mu\nu}=\partial_\mu\partial_\nu - \eta_{\mu\nu} \square$ that appears in the action, meaning that using the standard formula for a Gaussian integral is illegal (since it depends on the inverse of $O$). However, if we stick to sources which are conserved, with $\partial_\mu J^\mu=0$, then we restrict to a subspace where we can invert $O$. One way to see this is that in other gauges, you get a term like $p_\mu p_\nu$ in the propagator, but if $J$ is conserved then the Fourier transform $\tilde{J}$ obeys $p_\mu \tilde{J}^\mu=0$ and those gauge-dependent terms don't contribute to $W[J]$.
$W[J]$ is related to the interaction energy between the sources $J$ by (for example, see Section 1.4 of Zee's textbook, especially near Eq 6)
\begin{equation}
W = - E T
\end{equation}
where $T$ is the time interval of the experiment. Therefore the sign of $W$ determines whether the interaction energy is attractive or repulsive. If the energy between two like charges is negative then the force will be attractive, and if it is positive then the force will be repulsive.
We can factor the $\eta_{\mu\nu}$ factor out of the integral. For two static point charges with the same charge $q$ at positions $\vec{x}_1$ and $\vec{x}_2$, we have
\begin{equation}
J^0(x) = q\left(\delta^3(\vec{x}-\vec{x}_1) + \delta^3(\vec{x}-\vec{x}_2)\right), \ \ J^i(x) = 0
\end{equation}
and the interaction energy becomes
\begin{equation}
E = \left(q^2 \eta_{00}\right) \times \left[ - \int \frac{d^3 p}{(2\pi)^3} \frac{e^{i \vec{p} \cdot (\vec{x}_1 - \vec{x}_2)}}{\vec{p}^2} \right]
\end{equation}
The term in brackets is negative -- again see Zee's book. Therefore, the overall energy is positive, because the term in square brackets is multiplied by $\eta_{00}=-1$. This indicates that the force between like charges is repulsive. We also see that for opposite sign forces, the last term would yield a minus sign (because $J^0$ for the two charges would have opposite signs), yielding to an attractive force.
An important conclusion from this analysis is that Lorentz invariance is crucial for determining whether the force is attractive or repulsive, because of the factor of $\eta_{\mu\nu}$. This factor is missing for scalar fields, which is why scalars give rise to attractive forces, and this factor effectively appears squared for gravity because of the spin-2 nature of the graviton, which is why gravity is universally attractive.
