Here is a proof that goes along the lines of the OP's argument. The method is for numeric matrices, that is $x_i,y_j\in\mathbb{C}$ since it uses simple Calculus methods.
Notice that $$c_n(x_1,\ldots,x_n,y_1,\ldots,y_n):=\operatorname{det}\big(\frac{1}{x_i+y_j}\big)$$
is homogenous of over $-n$, that is
$$c_n(\lambda x_1,\ldots,\lambda x_n,\lambda y_1,\ldots,\lambda y_n)=\lambda^{-n}c_n(x_1,\ldots, x_n,y_1,\ldots,y_n)$$
For simplicity, set $\mathbf{x}=(x_1,\ldots, x_n)$ (similarly for $\mathbf{y}$). From the definition of determinant it follows that
$$
c_n(\mathbf{x},\mathbf{y})=\sum_{\sigma\in S_n}(-1)^{\sigma}\frac{1}{x_1+y_{\sigma(1)}}\cdot\ldots\cdot\frac{1}{x_n+y_{\sigma(n)}}=\frac{P(\mathbf{x},\mathbf{y})}{\prod_{1\leq i,j\leq n}(x_i+y_i)}
$$
where $P$ is a polynomial on $\mathbf{x}$ and $\mathbf{y}$ which is homogeneous of order $n^2-n$ (observe that $\prod_{1\leq i,j\leq n}(x_i+y_i)$ is the common denominator of all the rational expressions in the sum that defines $c$). Notice that if $x_\ell=x_k$ (or $y_\ell=y_k$) for some $\ell<k$, then the matrix $C_n=\big(\frac{1}{x_i+y_j}\big)$ would have rows (resp. columns) $\ell$ and $k$ identical. It follows that
$$c_n(\mathbf{x},\mathbf{y})=k_n\frac{\prod_{1\leq i<j\leq n}(x_i-x_j)(y_i-y_j)}{\prod_{1\leq i,j\leq n}(x_i+y_j)}$$
for some constant $k_n$.
I tried to find a particular choice of $\mathbf{x}$ and $\mathbf{y}$ that yields a matrix which a determinant easy to compute and which allow us to estimate $k_n$, but did not go very far. Then I turn to some simple Calculus: the function $\mathbf{x}\mapsto x_1 c(\mathbf{x},\mathbf{y})$ is the determinant of the matrix obtained from $C$ by multiplying the first row of $C$ by $x_1$ and leaving the other ones the same. The continuity of the determinant function yields
\begin{align}
\lim_{y_1\rightarrow\infty}\Big(\lim_{x_1\rightarrow\infty}x_1c_n(\mathbf{x},\mathbf{y})\Big)&=c_{n-1}(x_2,\ldots,x_n,y_2,\ldots,y_n)\\
&=k_{n-1}\frac{\prod_{2\leq i<j\leq n}(x_i-x_j)(y_i-y_j)}{\prod_{2\leq i,j\leq n}(x_i+y_j)}
\end{align}
On the other hand,
\begin{align}
\lim_{y\rightarrow\infty}\Big(\lim_{x_1\rightarrow\infty} x_1c_n(\mathbf{x},\mathbf{y})\Big)&=\lim_{y_1\rightarrow\infty}\Big(\lim_{x\rightarrow\infty}k_n\frac{x_1\prod_{1\leq i<j\leq n}(x_i-x_j)(y_i-y_j)}{\prod_{1\leq i,j\leq n}(x_i+y_j)}\Big)\\
&=k_n\frac{\prod_{2\leq i<j\leq n}(x_i-x_j)(y_i-y_j)}{\prod_{2\leq i,j\leq n}(x_i+y_j)}
\end{align}
Putting things together, we have that $k_n=k_{n-1}$
Clearly $k_1=1$. The desired formula follows.
Best Answer
This doesn't really have anything to do with Jacobians, it holds for any matrix.
If $j=k$, it's just expansion of the determinant along the $j$th row. And if $j \neq k$, it's “expansion along the wrong row”, where you expand along one row but use the cofactors from another row, which gives the same result as if you do expansion of a determinant where two rows are equal, i.e., the result is zero.