Write $S_k = \sum_{n \in \mathbb{Z}} (z-n)^{-k}$, summations taken symmetrically in the Eisenstein sense. By uniform/normal convergence, we know that $S_k' = -kS_{k+1}$. We are interested in $S_aS_b - S_c$.
By various theorems on uniform/normal convergence $S_k$ has a Laurent expansion with principal part $z^{-k} $ at $z=0$. Consideration of the expansion of $S_1$ (namely $S_1 =\pi\cot{\pi z} = \frac{1}{z} + \sum_{m=1}^{\infty} -2\zeta(2m) z^{2m-1}$) implies that $S_{2k} - z^{2k} \sim 2\zeta(2k)$ for $a=2k$ even and $-(2k-1)\zeta(2k)z$ for $a=2k-1$ odd.
It follows immediately that $S_a S_b \sim z^{-a-b} $ as $z \to 0$, so certainly we must have $a+b=c$. Are there any others than $a=1,b=2$ and $c=a+b=3$? First check $S_a^2$. The Laurent expansion begins
$$ z^{-2a}+2\alpha z^{-a}+\dotsb $$
or
$$ z^{-2a}+2\alpha z^{1-a}+\dotsb. $$
$a \geq 1$, and for $a=1$, we know that $S_1^2 = S_2-\pi^2$, so it suffices to check $a>1$. But then $-a$ and $1-a$ are both negative, so there is no way to cancel these out.
So now suppose that $a>b \geq 1$. The beginning of the principal part of the expansion of $S_aS_b$ is given by one of the four cases
$$ (z^{-a}+\alpha+O(z))(z^{-b}+\beta+O(z)) = z^{-a-b} + \beta z^{-a} + o(z^{a}) \\
(z^{-a}+\alpha z+O(z^2))(z^{-b}+\beta+O(z)) = z^{-a-b} + \beta z^{-a} + o(z^{a}) \\
(z^{-a}+\alpha+O(z))(z^{-b}+\beta z+O(z^2)) = z^{-a-b} + \beta z^{1-a} + \alpha z^{b} + o(z^{a}) \\
(z^{-a}+\alpha z+O(z^2))(z^{-b}+\beta z+O(z^2)) = z^{-a-b} + \beta z^{1-a} + o(z^{1-a})
$$
depending on the parity of $a$ and $b$. The first two of these are trivially disposed of. $a>1$, so we can also forget about the fourth. The third requires $a=b+1$ to be nontrivial. But $S_a S_{a+1}$ is the derivative of $(S_a)^2/(2a)$, so it suffices to check that we have no relation of the form $S_a^2 - S_{2a}+C$ except for $S_1^2 = S_2 - \pi^2$. But this again follows from the principal parts: there is a $2\alpha z^{-a}$ term or $2\alpha z^{1-a}$ that is not cancelled by the $S_{2a}$.
Hence the only possible identity of the given form is $S_1S_2=S_3$.
This is a general method that can be employed: we have the result
If $f$ is meromorphic and periodic with real period $L$, with finitely many poles $a_i$ in one period strip, and $f(z) $ converges to finite, but possibly different numbers as $\Im(z) \to \pm \infty$, then there are constants $c_{i,k}$ and $A$ so that
$$f(z) = A+\sum_{i} \sum_{k=1}^{n_i} c_{i,k}\cot^k{\frac{z-a_i}{L}}$$
(Whittaker and Watson has a weaker version for simple poles, but the proof can be easily effected by making the principal parts in the sum match those in $f$'s Laurent series. The different is then a bounded analytic function, which is constant by Liouville.)
Differentiating $S_1^2+\pi^2=S_2$ gives $-2S_1S_2 = -2S_3$, which is your relation. Differentiating again gives
$$3S_4 = S_2^2+S_1S_3 = S_2^2+S_1^2 S_2 = \pi^2+4\pi^2S_1^2+3S_1^4, $$
and so on. This type of argument, together with an induction, shows that any $S_k$ can be expressed as a polynomial in $S_1$ of degree $k$. This also means that there is an expansion of $f$ from the paragraph above of the form
$$ f(z) = A' +\sum_{i} \sum_{k=1}^{n_i} c_{i,k}' S_k\left(\frac{z-a_i}{L}\right).$$
This was all worked out in detail by Eisenstein in the nineteenth century, but I don't have the references to hand. A modern account is given in Walker's Elliptic funtions: A constructive approach, which is astonishingly hard to get hold of considering that it's only 20 years old.
