Let $$g(y)= \int_{\mathbb{R}^3}\frac{f(x)}{|x−y|}dx$$
Using Fourier theory or the argument here: Can convolution of two radially symmetric function be radially symmetric? $g$ will be spherically symmetric.
So, it is only necessary to find the value of $g$ for $y$ along the north pole: $y = (0,0,r)$; other $y$ with same magnitude $r$ will have the same value by spherical symmetry.
Setting $x=(s\sin(θ)\cos(ϕ),s\sin(θ)\sin(ϕ),s\cos(θ))$, observe that $|x-y| =\sqrt{r^2+s^2-2rs\cos(\theta)}$ for $y$ along the north pole.
So,
$$g(r)=\int_0^\infty s^2\int_0^{\pi}\sin(\theta)\int_0^{2\pi}\frac{f(s)}{\sqrt{r^2+s^2-2rs\cos(\theta)}}d\phi d\theta ds$$
$$\implies g(r)=2\pi\int_0^\infty f(s)s^2\int_0^{\pi}\sin(\theta)\frac{1}{\sqrt{r^2+s^2-2rs\cos(\theta)}} d\theta ds$$
$$\implies g(r)=2\pi\int_0^\infty f(s)\frac{r}{s}(r+s-|r-s|)ds$$.
So, we have reduced the 3-D convolution to a 1-D integral operator with kernel $K(r,s)=2\pi\frac{r}{s}(r+s-|r-s|)$. Obviously, there's no further simplification without knowing the form of $f$.