One needs to be careful when taking the $T\rightarrow 0$ limit. Consider $\epsilon_k = 0$. In this case, the integral form diverges. Therefore, the statement that this integral is always $0$ is clearly incorrect.
One takes a clue of the proper procedure by realizing that $\delta(\epsilon_k)$ is a distribution, and hence must be "properly" understood within an integral. Therefore, we really want to be considering the quantity:
$$\lim_{\beta\rightarrow\infty} \int d\epsilon_k f(\epsilon_k) \frac{1}{\beta}\sum_{n}\frac{1}{(i\omega_n - \epsilon_k)^2}$$
where we will take $f(\epsilon_k)$ to be as nice as possible (smooth, vanishes sufficiently quickly at $\infty$, etc). Now, we need to be painstakingly careful about the order of operations.
First, let's evaluate this quantity in the order shown.
\begin{align}
&\lim_{\beta\rightarrow\infty} \int d\epsilon_k f(\epsilon_k) \frac{1}{\beta}\sum_{n}\frac{1}{(i\omega_n - \epsilon_k)^2} \\
=& \lim_{\beta\rightarrow\infty} \int d\epsilon_k f(\epsilon_k) n_F'(\epsilon_k) \\
=& f(0)
\end{align}
Instead, we could have "moved the limit inside the integral" by treating $n_F(\epsilon_k)$ as a distribution, so that $\lim_{\beta\rightarrow \infty} n_F'(\epsilon_k) = \delta(\epsilon_k)$ in the sense of distributions, and arrived at the same result.
Now, how should we understand the limiting procedure converting the sum into an integral in the sense of distributions? Using our "sufficiently nice" properties of $f(\epsilon_k)$, we can interchange the sum and the integral:
\begin{align}
&\lim_{\beta\rightarrow\infty} \frac{1}{\beta}\sum_{n} \int d\epsilon_k f(\epsilon_k) \frac{1}{(i\omega_n - \epsilon_k)^2} \\
=& \int \frac{d\omega}{2\pi} \int d\epsilon_k f(\epsilon_k) \frac{1}{(i\omega - \epsilon_k)^2}
\end{align}
where in the second line we used the definition of the (Riemann) integral. We wish to interchange these two integrals so that we evaluate $\int d\omega$ first; classically this is not allowed since the integral is not absolutely convergent. However, let us do the following trick:
\begin{align}
&\int \frac{d\omega}{2\pi} \int d\epsilon_k f(\epsilon_k) \frac{1}{(i\omega - \epsilon_k)^2} \\
=& \lim_{\epsilon \rightarrow 0} \int_{\mathbb{R}-(-\epsilon,\epsilon)} \frac{d\omega}{2\pi} \int d\epsilon_k f(\epsilon_k) \frac{1}{(i\omega - \epsilon_k)^2} \\
=& \lim_{\epsilon \rightarrow 0} \int d\epsilon_k f(\epsilon_k) \int_{\mathbb{R}-(-\epsilon,\epsilon)} \frac{d\omega}{2\pi} \frac{1}{(i\omega - \epsilon_k)^2} \\
=& \lim_{\epsilon \rightarrow 0} \int d\epsilon_k f(\epsilon_k) \frac{1}{\pi}\frac{\epsilon}{\epsilon^2 + \epsilon_k^2} \\
=& f(0)
\end{align}
In the second line, we remove the problematic region so that we can interchange the integrals; we take the principal value of the integral. Note that in the second to last line, a nascent delta-function $\frac{1}{\pi}\frac{\epsilon}{\epsilon^2 + \epsilon_k^2}$ appears!
In other words, in order to take the $\beta\rightarrow \infty$ limit first to convert the sum into an integral, we need to be calculating the principal value of the integral, with the $\epsilon \rightarrow 0$ limit being taken in the sense of distributions (ie after the $\int d\epsilon_k$ integral). Concretely, this means
\begin{align}
&\frac{1}{\beta}\sum_n\frac{1}{(i\omega_n -\epsilon_k)^2}\\
=&\int_{\mathbb{R}-(-\epsilon,\epsilon)} \frac{d\omega}{2\pi} \frac{1}{(i\omega - \epsilon_k)^2} \\
=& \int_{\mathbb{R}-(-\epsilon,\epsilon)} \frac{d\omega}{2\pi} \frac{1}{(i\omega - \epsilon_k)^2} \\
=& \frac{1}{\pi}\frac{\epsilon}{\epsilon^2 + \epsilon_k^2} \\
=& \delta(\epsilon_k)
\end{align}
