Which limit for Matsubara frequency summation?
Yes, for this particular Matsubara sum you mentioned, taking different limits will lead to two results differed by 1. This is because the summation in consideration does not converge, which can be seen from the following integral (the continuous limit of the sum) considering the ultraviolet divergence (the large $\omega$ behavior) $$\frac{1}{2\pi}\int_{-\infty}^{\infty}\mathrm{d}\omega\frac{1}{i\omega-\epsilon}\sim\frac{1}{2\pi i}\int_{-\infty}^{\infty}\mathrm{d}\omega\frac{1}{\omega}\rightarrow\infty.$$ Asking for the result of an essentially divergent sum will not even lead to a definite answer. Introducing the factor $e^{i\omega_n\tau}$ is to control the convergence of the sum, but the result will then depend on how one chooses to control the convergence, i.e. $\tau\rightarrow0^+$ (controlling the left-half complex plane) or $0^-$ (controlling the right-half complex plane).
The two results will be different by a shift of 1. Even this difference "1" has a physical meaning, reflecting the ambiguity in defining Fermion particle number. Because the summation $$\frac{1}{\beta}\sum_n\frac{1}{i\omega_n-\epsilon}=G(\tau=0)=-\mathcal{T}_{\tau}\langle\psi(\tau=0)\bar{\psi}\rangle$$ corresponds to the Green's function at imaginary time 0, which physically means to count the number of Fermions on the energy level $\epsilon$. Here $\mathcal{T}_{\tau}$ stands for the time ordering operator. However $\tau=0$ can be either understood as $\tau\rightarrow0^+$ or $\tau\rightarrow0^-$ (corresponding to the two ways to control the convergence), which becomes very critical here because the time ordering operator will then order the two cases differently: $$-\mathcal{T}_{\tau}\langle\psi(\tau=0^+)\bar{\psi}\rangle = -\langle\psi\bar{\psi}\rangle,$$ $$-\mathcal{T}_{\tau}\langle\psi(\tau=0^-)\bar{\psi}\rangle = \langle\bar{\psi}\psi\rangle.$$ According to the anticommutation relation of the Fermion operators, we have $\bar{\psi}\psi=1-\psi\bar{\psi}$, that is why it is expected that the two results will differ by a constant 1. In fact $\bar{\psi}\psi$ counts the number of particles, while $\psi\bar{\psi}$ counts the number of holes. It becomes a problem whether to define the Fermion number (with respect to the vacuum state) as the particle number or as the negative of the hole number. This depends on our definition of the vacuum state: whether we treat the vacuum as no particle state or as a state filled with particles (like the concept of Fermi sea). We know both choices are acceptable, just as we can either define electron as a particle in the empty space or as a hole in the anti-electron Fermi sea. Because of this ambiguity, one needs to specify the choice when the $\tau\rightarrow0$ limit is taken. In conclusion, both ways of taking the limit are legitimate, and the difference in the results is just a matter of the definition of Fermion vacuum state.
Similar discrepancy also exists in the Bosonic case, which is again related to the definition of Boson vacuum state. In general, such discrepancy roots in the divergence of the Matsubara summation in consideration. If the summation itself converges to a definite result, then the result will be unique no matter which limit is chosen to take. For example, if one tries to calculate $$\frac{1}{\beta}\sum_n\frac{1}{(i\omega_n-\epsilon)^2},$$ the result will always be $\beta/4 \mathrm{sech}^2 (\beta\epsilon/2)$ (for Fermionic case) no matter $\tau\rightarrow0^+$ or $0^-$.