You are absolutely right that the limit in which this approximation holds is
$$\beta(\epsilon - \mu) \gg 1 \,,$$
which is not trivially the 'high-temperature limit', and indeed looks rather like the low temperature limit. However, it also looks like the limit of large negative $\mu$. If we want to know how temperature will affect the exponent, we need to know how temperature will effect the chemical potential. To proceed, suppose we're dealing with a gas of non-interacting particles. The grand potential is, in this limit,
$$ \Phi = -k_B T \ln \mathcal{Z} = -k_B T \int_0^\infty \ln \mathcal{Z}_\epsilon \,g(\epsilon)\,\mathrm{d}\epsilon \simeq -k_B T \int_0^\infty \ln \bigg(1 + \exp(-\beta(\epsilon - \mu))\bigg)\,g(\epsilon) \,\mathrm{d}\epsilon \,,$$
where $\mathcal{Z}_\epsilon$ is the grand partition function associated with the energy level $\epsilon$ and $g(\epsilon)$ is the density of states. The integral is essentially just a sum of the partition functions due to each energy level. To get to the final expression we have assumed that we can approximate the grand partition function like so:
$$ \mathcal{Z}_\epsilon = \sum_{n} \bigg(\exp(-\beta(\epsilon - \mu))\bigg)^n \simeq 1 + \exp(-\beta(\epsilon - \mu)) \,,$$
which corresponds to the limit stated at the top. As a brief detour, if we want to find the average occupancy of the energy level $\epsilon$, we can use
$$ \langle N_\epsilon \rangle = -\left(\frac{\partial \Phi_\epsilon}{\partial \mu}\right)_{T,V} \simeq \exp(-\beta(\epsilon - \mu))\qquad \mathrm{where} \qquad \Phi_\epsilon = -k_B T \ln \mathcal{Z}_\epsilon\,,$$
which is the Maxwell-Boltzmann distribution we were expecting (in the second equality we have Taylor expanded the logarithm in accordance with $\beta(\epsilon - \mu) \gg 1$). Now the density of states for a three-dimensional gas in a box can be obtained by standard means --- I won't bother going through it here, but the end result is:
$$ \Phi = -k_B TV\left(\frac{mk_B T}{2 \pi \hbar^2}\right)^{3/2} \exp(\beta \mu) \equiv -\frac{k_B T V}{\lambda^3} \exp(\beta \mu) \,,$$
where the thermal wavelength $\lambda$ has been defined appropriately. From here we can write
$$N_\mathrm{tot} \equiv N = -\left(\frac{\partial \Phi}{\partial \mu}\right)_{T,V} = \frac{V}{\lambda^3} \exp(\beta \mu) \,,$$
and hence
$$ \boxed{\mu = k_B T \ln \left(\frac{N \lambda^3}{V}\right) \,.}$$
Now to answer your question. The condition at the top can be considered the limit of $\beta \mu$ being large and negative. We see from the above that
$$ \beta \mu = \ln\left(\frac{N \lambda^3}{V}\right) \qquad \mathrm{where} \qquad \lambda = \left( \frac{2 \pi \hbar^2}{mk_B T}\right)^{1/2} \,.$$
This quantity will be large and negative when the argument of the logarithm is small. This will be the case for a) low densities $N/V$, b) high temperatures $T$ and/or c) high-mass particles.
You should think of the underlying situation in which the classical limit holds as when the number of thermally accessible states vastly exceeds the number of particles. This is because under such circumstances we can ignore multiple occupation of energy levels, which means we can ignore the fine details of particle indistinguishability. In the canonical distribution, when the number of states vastly exceeds the number of particles we can account for indistinguishability with a simple (but approximate) correction of division of the partition function by $N!$ --- we must do this even in the classical case, otherwise we run into all sorts of problems like the Gibbs paradox. However, when states start to become multiply occupied, this simple prescription fails, and we need to be more sophisticated in our consideration of particle indistinguishability.
If you imagine our gas particles as being wavepackets with a width of $\lambda$ as defined above, then you can think of each particle as occupying a volume $\lambda^3$. This has a nice interpretation --- the quantity $N \lambda^3/V$ that appears in the expression for the chemical potential can be thought of as the fraction of space occupied by the particles. The classical limit corresponds to this quantity being small, so that it's very unlikely for two particles to be in the same place --- i.e., be in the same state (here I'm essentially considering the states of our system to be position eigenstates rather than the usual energy eigenstates). If this quantity becomes larger, we start to get 'multiple-occupation', and so we imagine our classical approximation will break down. This is consistent: when $N \lambda^3 /V \sim 1$, the argument of the logarithm in the chemical potential is no longer large and negative, and so indeed the condition at the very top of this page breaks down.
Hope this helps!
Best Answer
It takes a lengthy proof, but Lyman Spitzer shows in the second chapter of Physical Processes in the Interstellar Medium (the standard text in interstellar matter studies) that the velocity distribution of interstellar gas particles (which is what forms nebulae) is very nearly Maxwellian - the deviation is less than 1%.
Other larger systems, probably not so much - Maxwell-Boltzmann statistics work best when kinetic energy is dominant in a system. But I don't know much about the topic, so that is a guess.