Elementary proof that $\sum_{j=1}^{n} \prod_{k \neq j} \frac{1}{1+(a_j - a_k + i)^2} \in \mathbb{R}$ for $a_1, \dots, a_n \in \mathbb{R}$ distinct

Since nobody seems to be willing to answer, let me expand my comment above.

Consider a more general problem. Suppose we started with $$\frac1{2\pi}\int_{-\infty}^{\infty}\frac{dx}{\prod_{k=1}^n(x-z_k)(x-\bar z_k)},$$ where all $z_k$'s are pairwise distinct and all $\Im z_k>0$. The analog of your problem would then be to show that $$\sum_{j=1}^n\frac{i}{(z_j-\bar z_j)\prod_{k\ne j}^n(z_j-z_k)(z_j-\bar z_k)}\; \in\mathbb R,$$ which is the same as to show that $$\sum_{j=1}^n\frac{1}{(z_j-\bar z_j)\prod\limits_{k\ne j}^n(z_j-z_k)(z_j-\bar z_k)}+\sum_{j=1}^n\frac{ 1}{(\bar z_j- z_j)\prod\limits_{k\ne j}^n(\bar z_j-z_k)(\bar z_j-\bar z_k)}=0\tag{1}$$ This is a special case of the identity $$\qquad\qquad\sum_{j=1}^N\frac{1}{\prod_{k\ne j}^N(u_j-u_k)}=0,\tag{2}$$ which corresponds to setting in (2) $N=2n$ and $u_j=z_j$, $u_{j+n}=\bar z_j$ for $j=1,\ldots,n$. It thus remains to prove (2).

Knowing the beginning of the story, one proof is straightforward: it suffices to compute $\displaystyle\int_\Gamma\frac{du}{\prod_{k=1}^N(u-u_j)}$ by residues in two different ways (here $\Gamma$ is a closed contour around $u_1,\ldots,u_N$). Shrinking the contour $\Gamma$ produces the left side of (2), while expanding it to infinity gives zero.

If instead we want an algebraic proof of (2), we may use Lagrange interpolation. Given $f_1,\ldots,f_N\in\mathbb C$, there exists a unique polynomial $Q(x)$ of degree $N-1$ such that $Q(u_k)=f_k$, and moreover $$Q(x)=\sum_{j=1}^Nf_j\prod_{k\ne j}^N\frac{x-u_k}{u_j-u_k}.$$ Choose $f_1=\ldots=f_N=1$ so that $Q(x)=1$. The previous formula then transforms into $$\sum_{j=1}^N\prod_{k\ne j}^N\frac{x-u_k}{u_j-u_k}=1.$$ Equating the coefficients of $x^{N-1}$ at both sides immediately yields (2).