The point is that to achieve the sum over Matsubara frequencies
we can use a contour integral
$$\oint_C g(z) f(z)$$
with the contour described in fig 1 here, so long as we choose an $f(z)$ with simple poles exactly at the Matsubara frequencies $\omega_n$. This determines $f(z)$ to be proportional to the Bose-Einstein distribution. In fact, if we think about the sum as computed in a correlation function, we could have derived the integral above instead by considering computing the expectation value in the thermal density matrix $\int dE n(\beta,E) |E\rangle\langle E|$ determined by the Bose-Einstein distribution $n(\beta,E)$. I think this is an equally good starting point of the logic, then deriving the Matsubara frequencies as the poles of the distribution function.
I'd say the whole point of the Matsubara frequencies is that they are the poles of the distribution function.