consider the grand canonical ensemble,
$$ \rho \sim \exp[-\beta (E-\mu N)] $$
In the exponent, the inverse temperature $\beta = 1/kT$ is the coefficient in front of one conserved quantity, the (minus) energy, while another coefficient, $\beta\mu$, is in front of the number of particles $N$. The chemical potential is therefore the coefficient in front of the number of particles - except for the extra $\beta$. The number of particles has to be conserved as well if $\mu$ is non-zero.
As long as the sign of $N$ is well-defined, the sign of the chemical potential is well-defined, too.
Now, for bosons, $\mu$ can't be positive because the distribution would be an exponentially increasing function of $N$. Note that in the grand canonical ensemble - which is really the ensemble in which $\mu$ is sharply defined - the dual variable to $\mu$, namely $N$, is not sharply defined. However, the probability for ever larger - infinite $N$ would be larger, so the distribution would be peaked at $N=\infty$. Such a distribution couldn't be well-defined. We want, in the thermodynamic limit, the grand canonical ensemble, while assuming a fixed $\mu$, also generate a finite and almost well-defined $N$, within an error margin that goes to zero in the thermodynamic limit. That couldn't happen for bosons and a positive $\mu$.
This catastrophe would be possible because $N=\sum_i n_i$, a summation over microstates $i$, and every $n_i$ can be an arbitrarily high integer for bosons. For fermions, the problem doesn't occur because $n_i=0$ or $1$ for each state $i$. So for fermions, we can't argue that $\mu$ has to be positive. Note that $E-\mu N$ in the exponent is the sum $\sum_i N_i (e_i-\mu)$. For bosons, the problem occurred for states for which $e_i-\mu$ was negative i.e. $e_i$ was low enough. For fermions, however, the number of such states - and therefore the maximum number of fermions in them - is finite, so the divergence doesn't occur if the chemical potential is positive. For fermions, positive $\mu$ is OK.
In fact, for fermions, both positive and negative $\mu$ is OK. Also, it is easy to see that if both particles and antiparticles exist, $\mu$ of the antiparticle has to be minus $\mu$ of the particle because only the difference $N_{\rm particles} - N_{\rm antiparticles}$ is conserved; this is true both for bosons and fermions.
So if the potential for electrons is positive, the potential for positrons or holes (which play the very same role) has to be negative, and vice versa.
Nothing changes about the meaning of the chemical potential when one switches from classical physics to quantum physics: in fact, above, I was assuming that there are "discrete states" for the particles, just like in quantum physics - otherwise we wouldn't be talking about bosons and fermions which are only relevant in the quantum setup. Classical physics is a limit of quantum physics in which the number of states is infinite because $\hbar$ goes to zero, so a finite number of particles never ends up in "exactly the same state". In some sense, Ludwig Boltzmann, while working in the context of classical statistical physics, was inherently using the thinking and intuition of quantum statistical physics - he was a truly ingenious "forefather" of quantum physics.
In relativity, one has to be careful how we define the energy of a state. Note that the physically meaningful combination that appears in the exponent is $e_i-\mu$, so if one shifts $e_i$ e.g. by $mc^2$, the latent energy, one has to shift $\mu$ in the same direction by the same amount. The notions of chemical potential obviously work in relativity, too. Relativistic physics is not a "completely new type of physics". It is just a type of the old physics that happens to respect a symmetry - the Lorentz symmetry.
Again, in quantum field theory which combines both quantum mechanics and relativity, statistical physics including the notion of the chemical potential also works but one must be careful that particle-antiparticle pairs may be created with enough energy. That implies $\mu_{p}=-\mu_{np}$, as I said.
There can't be any general ban on a negative chemical potential of fermions: fermionic $\mu$ can have both signs. However, in the particular theory that Allan wanted to describe, he could have had more detailed reasons why $\mu$ should have been positive for his fermions. I am afraid that this would be an entirely different, more specific question - one about superconductors. As stated, your question above was one about statistical physics and I tend to believe that the text above exhausts all universal facts about the sign of the chemical potential in statistical physics.
Throughout, let's assume that the ground state energy of the system under consideration is zero.
Chay Paterson has addressed your question in the case of a gas of bosons in which the number of particles is not conserved, but from the wording of your question, it seems that you're concerned about the case in which the total number of particles is fixed.
For a system of bosons with a fixed number $N$ of particles, the answer to your question is
No. The chemical potential is nonzero for all $T>0$.
You point out, however, that it is often said that below the critical temperature $T_c$, the chemical potential is zero, so what's going on? The resolution is essentially that
the chemical potential very well-approximated by zero for almost all temperatures below the critical temperature, but it is never exactly zero.
As you cool the system down from above the critical temperature, the number $N_e$ of particles in excited states get's lower and lower. On the other hand, the number of particles in the excited states at any given temperature is bounded above, and this bound decreases as a function of temperature like $T^{3/2}$ (for non-relativistic systems in three dimensions). At a particular sufficiently low temperature (the critical temperature), this upper bound is lower than the total number of particles in the system and particles are forced into the ground state. But at this point, the chemical potential does not drop exactly to zero. It does, however, get very small very quickly as the temperature decreases and the number of particles in the ground state increases.
In fact, If we look at the number of particles in the ground state as a function of temperature
$$
N_0 = \frac{1}{e^{-\mu(T)/kT}-1}
$$
which gives
$$
\mu(T) = -kT\ln\left(1+\frac{1}{N_0}\right)
$$
then we see that the chemical potential is always strictly less than zero because the argument of the log is always strictly greater than $1$. But as you dial the temperature down below the critical temperature, and as the number of particles in the ground state increases towards $N$, the total number of particles in the system, which is presumably very large, and the chemical potential decreases since the argument of the log approaches $1$.
Best Answer
To determine the upper limit on chemical potential for a gas of $\mathcal N$ bosons, look at the form of the Bose distribution in the grand canonical ensemble with $\langle N \rangle = \mathcal N$. When using the GCE, it's easiest to work at chemical potential $\mu$ and to then choose $\mu(\mathcal N)$ so that $\langle N\rangle(\mu)=\mathcal N$. Each state $s$ has average occupancy $$ \langle n_s\rangle=\frac{\sum_{n\geq 0} ne^{-\beta n(\epsilon_s-\mu)}}{\sum_{n\geq 0} e^{-\beta n(\epsilon_s-\mu)}}=\frac{1}{\Xi_s}\frac{\partial}{\partial(\beta\mu)}\Xi_s,\quad \Xi_s=\frac{1}{1-e^{-\beta(\epsilon_s-\mu)}},\\ =-\partial_{(\beta\mu)}\log(1-e^{-\beta\epsilon_s+(\beta\mu)})=\frac{e^{\beta\mu}}{1-e^{-\beta(\epsilon_s-\mu)}}. $$ This is finite as long as $\mu<\epsilon_s$. In order for $\langle N\rangle=\sum_s \langle n_s\rangle$ to be finite we need $\mu<\min_s \epsilon_s=\epsilon_0$. Hence, for any system of bosons where $N$ is conserved we have $\mu<\epsilon_0$. It is conventional to set $\epsilon_0=0$ for simplicity, but you can have systems with $\epsilon_0\neq 0$. As you implied by your final question, in these systems the critical value of $\mu$ is $\epsilon_0$ in the thermodynamic limit, with $N, V, E\rightarrow \infty$ and $\mu, p, T$ held constant. Of course, if the system has no BEC phase then as $T\rightarrow 0$, the chemical potential $\mu$ never exceeds some value $\mu_\max<\epsilon_0$.