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^-$.
I) The Euclidean path integral with $\hbar=1$ reads
$$ Z~=~\int_{DBC} \!{\cal D}x ~e^{-S},\tag{1}$$
with Dirichlet boundary conditions (DBC)
$$ x(0)~=~0~=~x(T).\tag{2}$$
We expand the real periodic variable $x\in\mathbb{R}$ in Fourier series$^1$
$$ \begin{align} x(t) ~=~&\frac{a_0}{2}+\sum_{n\in\mathbb{N}} \left\{a_n \cos(\omega_n t) + b_n \sin(\omega_n t)\right\}\cr
~=~& \sum_{n\in\mathbb{Z}} c_n e^{i\omega_n t}, \cr
\omega_n~:=~&\frac{2\pi n}{T},\qquad a_n,b_n~\in~ \mathbb{R},\cr
c_n~\in~& \mathbb{C}, \qquad c_n^{\ast}~=~c_{-n}. \end{align} \tag{3} $$
The DBC (2) becomes
$$ \sum_{n\in\mathbb{Z}} c_n~=~0\qquad\Leftrightarrow\qquad c_0 ~=~ -2{\rm Re}\sum_{n\in\mathbb{N}}c_n .\tag{4}$$
The action for a free non-relativistic point particle with mass $m=1$ reads:
$$ S~=~\frac{1}{2} \int_0^T \!dt~ \dot{x}^2~=~T\sum_{n\in\mathbb{N}}\omega_n^2 |c_n|^2 . \tag{5} $$
II) We know that the proper normalization of the path integral (1) is
$$ Z~=~\frac{1}{\sqrt{2 \pi T}}.\tag{6} $$
This can e.g. be deduced (without introducing fudge factors!) from the (semi)group property of Feynman path integrals, cf. this Phys.SE post and links therein. Up till now we have basically just restated what OP wrote in his question.
III) Now we would like to repeat the same calculation using Fourier series, i.e. work with the Matsubara frequencies. In this answer, we will not explore the (semi)group property, but just do a quick and dirty calculation using various fudge factors, and see what we get. Since this is homework, the explanation will be somewhat brief.
To make heuristic sense of the path integral (1), we will use the following zeta function regularization rules:
$$ \prod_{n\in\mathbb{N}} a ~=~\frac{1}{\sqrt{a}} \quad\text{and}\quad \prod_{n\in\mathbb{N}} n ~=~\sqrt{2\pi}, \tag{7}$$
stemming from the zeta function values
$$ \zeta(0)~=~-\frac{1}{2}
\quad\text{and}\quad
\zeta^{\prime}(0)~=~-\ln\sqrt{2\pi} , \tag{8} $$
respectively. Now let the path integral measure be
$$ \begin{align} {\cal D}x~:=~&\delta\left(B\sum_{n\in\mathbb{Z}} c_n\right) A\mathrm{d} c_0 \prod_{n\in\mathbb{N}} A^2 \mathrm{d}^2c_n \cr
~\stackrel{(7)}{=}~&\frac{1}{B}\delta\left(\sum_{n\in\mathbb{Z}} c_n\right) \mathrm{d} c_0 \prod_{n\in\mathbb{N}} \mathrm{d}^2c_n , \end{align}\tag{9} $$
where $A$, $B$ are fudge factors. Interestingly, eq. (9) is independent of the $A$-fudge factor! After performing the delta function integration and the Gaussian integrals, we find
$$ \begin{align} Z~=~&\frac{1}{B} \prod_{n\in\mathbb{N}} \frac{\pi}{T\omega_n^2}\cr
~=~&\frac{1}{B} \prod_{n\in\mathbb{N}} \frac{T}{4\pi n^2}\cr
~\stackrel{(7)}{=}~&\frac{1}{B\sqrt{\pi T}}. \end{align}\tag{10}$$
Apparently we should chose the fudge factor $B=\sqrt{2}$ in order to achieve the correct normalization (6).
--
$^1$ Note that the sine (cosine) modes (3) correspond trivially (non-trivially) to the even (odd) modes in my Phys.SE answer here, respectively.
Best Answer
From wikipedia, look at this contour:
This contour is the sum of a large circle at infinity and a tight contour going around the poles on the imaginary axis. The integration of the tight contour gives you the terms of the sum your interested in.
Along the large circle this term goes to 0:
$$ \frac{1}{z - (\epsilon - \mu)/\hbar} $$
So adding this contour to the tight one (the one that gives you the terms of the sum) is just adding zero. But now if you add the two contour together you get a simpler contour that just goes around the poles of the function your summing over. (The one depicted in this picture)
Therefore this equation:
$$\lim_ {\eta \downarrow 0 } \frac{1}{2 \pi i} \int _C dz \frac{e^{\eta z}}{z - (\epsilon - \mu)/\hbar} \frac{\mp 1}{e^{\hbar \beta z}\mp 1} = \lim_ {\eta \downarrow 0 }\frac{1}{\hbar \beta}\sum_{n}\frac{e^{i \omega_n \eta}}{i \omega_n - (\epsilon - \mu)/\hbar}.$$
still holds while using this contour:
Applying the residue theorem will then give you the contribution from the pole:
$$ \frac{1}{z - (\epsilon - \mu)/\hbar} $$ And thus you get: $$ \lim_ {\eta \downarrow 0 } e^{\eta z} \frac{\mp 1}{e^{\hbar \beta z}\mp 1}|_{z = (\epsilon - \mu)/\hbar}$$
Which is the result your looking for.