Fermi-Dirac and Bose-Einstein energy occupation number $n(\epsilon)$ in natural units ($[T]=[\epsilon]$) read $$n(\epsilon) = \frac{D(\epsilon)}{e^{(\epsilon-\mu)/T}\pm 1},$$ where $D(\epsilon)$ is the density of states at the given energy $\epsilon$, $+$ is Fermi-Dirac and $-$ is Bose-Einstein.
The usual answer to the question by which limit to arrive to Maxwell-Boltzmann statistics is: $$\frac{\epsilon-\mu}{T} \gg 1$$ or $$\frac{\epsilon_{min}-\mu}{T} \gg 1$$ This gives you purely formally Maxwell-Boltzmann distribution but is a make-it-work pseudo-argument as when $T \to \infty$, $\frac{\epsilon-\mu}{T} \to 0$ and on the contrary $T \to 0^+$, $\frac{\epsilon-\mu}{T} \to \infty$. Hence this rule would tell us to apply Boltzmann statistics to low temperatures and stick with the quantum statistics at high temperatures. This is obviously a fall into the fiery pit of slurry textbook mumbo jumbo.
To resolve this, I believe we have to presume a growing density of states $D(\epsilon)$ and faithfully take the $T \to \infty$ limit showing that every macroscopic feature (i.e. every moment of the original distribution) is in the limit reproduced by Maxwell-Boltzmann up to $\mathcal O(T^{-2},\sqrt{N},...)$.
The problem in the "textbook" argument as well as with the discussed limit is that $\epsilon-\mu$ actually always passes through values which are both greater and smaller than $T$ in the integration. The growing $T$ just "smears out" the distribution out into larger and larger regions in $\epsilon$, which motivated my guess that $D(\epsilon)$ must be growing so that the "smearing out" makes the high-energetic states dominate.
So how is the limit rigorously done? And are there some extra assumptions which aren't usually mentioned? (My guess is that the junk with $\mu \approx \epsilon$ in Bose-Einstein also needs some handling.)